Utility metric for unsupervised feature selection

Manifold learning considering non-linearities

Given is a data matrix $\mathbf{X}\in {\mathbb{R}}^{N\times d}$, with $\mathbf{X}=[{\mathbf{x}}_1;{\mathbf{x}}_2;\dots ;{\mathbf{x}}_N]$, ${\mathbf{x}}_i=[x_i^{(1)},x_i^{(2)},\dots ,x_i^{(d)}]$, $i=1,\dots ,N$, N the number of data points, and d the number of features (i.e., dimensions) in the data. The aim is to learn the structure hidden in the d-dimensional data and approximate it with only a subset of the original features. In this paper, this structure will be identified by means of clustering, where the dataset is assumed to be characterized by c clusters.

In spectral clustering, the clustering structure of this data can be obtained by studying the eigenvectors derived from a Laplacian built from the original data (Von Luxburg (2007), Biggs, Biggs & Norman (1993)). The data is represented using a graph $G=(\mathcal{V},\mathcal{E})$. $\mathcal{V}$ is the set of vertices v_i, i = 1,…,N where v_i = x_i. $\mathcal{E}=\{e_{ij}\}$ with i = 1,…,N j = 1,…,N is the set of edges between the vertices where {e_ij} denotes the edge between vertices v_i and v_j. The weight of these edges is determined by the entries w_ij ≥ 0 of a similarity matrix W. We define the graph as undirected. Therefore, the similarity matrix W, is symmetric (since w_ij = w_ji, with the diagonal set to w_ii = 0).

Typically, W is computed after coding the pairwise distances between all N data points. There are several ways of doing this, such as calculating the k-nearest neighbours (KNN) for each point, or choosing the ε-neighbors below a certain distance (Belkin & Niyogi, 2002).

In this paper, two similarity matrices are adopted inspired by the work in (Cai, Zhang & He, 2010), namely a binary one and one based on an RBF kernel. The binary weighting is based on KNN, being w_ij = 1 if and only if vertex i is within the K closest points to vertex j. Being a non-parametric approach, the binary embedding allows to simply characterize the connectivity of the data.

Additionally, the use of the RBF kernel is considered, which is well suited for non-linearities and allows to characterize complex and sparse structures (Von Luxburg, 2007). The RBF kernel is defined as K(x_i,x_j) = exp(−||x_i − x_j||²/2σ²). The selection of the kernel parameter σ is a long-standing challenge in machine learning. For instance, in Cai, Zhang & He (2010), σ² is defined as the mean of all the distances between the data points. Alternatively, a rule of thumb, uses the sum of the standard deviations of the data along each dimension (Varon, Alzate & Suykens, 2015). However, the estimation of this parameter is highly influenced by the amount of features or dimensions in the data, making it less robust to noise and irrelevant features. In the next section, a new and better informed method to approximate the kernel parameter is proposed.

The graph G, defined by the similarity matrix W, can be partitioned into multiple disjoint sets. Given the focus on multi-cluster data of our approach, the k-Way Normalized Cut (NCut) Relaxation is used, as proposed in Ng, Jordan & Weiss (2002). In order to obtain this partition, the degree matrix D of W must be calculated. D is a diagonal matrix for which each element on the diagonal is calculated as $D_{ii}={\sum }_j{W}_{i,j}$. The normalized Laplacian L is then obtained as L = D^−1/2WD^−1/2, as suggested in Von Luxburg (2007). The vectors y embedding the data in L can be extracted from the eigenvalue problem (Chung & Graham, 1997):

(1) $$\mathbf{L}\mathbf{y}=\lambda \mathbf{y}$$

Given the use of a normalized Laplacian for the data embedding, the vectors y must be adjusted using the degree matrix D:

(2) $$\alpha ={\mathbf{D}}^{1/2}\mathbf{y},$$

which means that α is the solution of the generalized eigenvalue problem of the pair W and D. These eigenvectors α are a new representation of the data, that gathers the most relevant information about the structures appearing in the high-dimensional space. The c eigenvectors, corresponding to the c highest eigenvalues (after excluding the largest one), can be used to characterize the data in a lower dimensional space (Ng, Jordan & Weiss, 2002). Thus, the matrix $\mathbf{E}=[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_c]$ containing column-wise the c selected eigenvectors, will be the low-dimensional representation of the data to be mimicked using a subset of the original features, as suggested in Cai, Zhang & He (2010).

Kernel parameter approximation for high-dimensional data

One of the most used similarity functions is the RBF kernel, which allows to explore non-linearities in the data. Nevertheless, the kernel parameter σ² must be selected correctly, to avoid overfitting or the allocation of all data points to the same cluster. This work proposes a new approach to approximate this kernel parameter, which will be denoted by ${\hat{\sigma }}^2$ when derived from our method. This method takes into account the curse of dimensionality and the potential irrelevant features or dimensions in the data.

As a rule of thumb, σ² is approximated as the sum of the standard deviation of the data along each dimension (Varon, Alzate & Suykens, 2015). This approximation grows with the number of features (i.e., dimensions) of the data, and thus, it is not able to capture its underlying structures in high-dimensional spaces. Nevertheless, this σ² is commonly used as an initialization value, around which a search is performed, considering some objective function (Alzate & Suykens, 2008; Varon, Alzate & Suykens, 2015).

The MCFS algorithm skips the search around an initialization of the σ² value by substituting the sum of the standard deviations by the mean of these (Cai, Zhang & He, 2010). By doing so, the value of σ² does not overly grow. This estimation of σ² suggested in Cai, Zhang & He (2010) will be referred to as ${\sigma }_0^2$. A drawback of this approximation in high-dimensional spaces is that it treats all dimensions as equally relevant for the final estimation of σ²₀, regardless of the amount of information that they actually contain.

The aim of the proposed approach is to provide a functional value of σ² that does not require any additional search, while being robust to high-dimensional data. Therefore, this work proposes an approximation technique based on two factors: the distances between the points, and the number of features or dimensions in the data.

The most commonly used distance metric is the euclidean distance. However, it is very sensitive to high-dimensional data, deriving unsubstantial distances when a high number of features is involved in the calculation (Aggarwal, Hinneburg & Keim, 2001). In this work, the use of the Manhattan or taxicab distance (Reynolds, 1980) is proposed, given its robustness when applied to high-dimensional data (Aggarwal, Hinneburg & Keim, 2001). For each feature l, the Manhattan distance δ_l is calculated as:

(3) $${\delta }_l={\displaystyle \frac{1}{N}\sum _{i,j=1}^N|x_{il}-x_{jl}|}$$

Additionally, in order to reduce the impact of irrelevant or redundant features, a system of weights is added to the approximation of ${\hat{\sigma }}^2$. The goal is to only take into account the distances associated to features that contain relevant information about the structure of the data. To calculate these weights, the probability density function (PDF) of each feature is compared with a Gaussian distribution. Higher weights are assigned to the features with less Gaussian behavior, i.e., those the PDF of which differs the most from a Gaussian distribution. By doing so, these will influence more the final ${\hat{\sigma }}^2$ value, since they allow a better separation of the structures present in the data.

Figure 1 shows a graphical representation of this estimation. The dataset in the example has 3 dimensions or features: f₁, f₂ and f₃. f₁ and f₂ contain the main clustering information, as it can be observed in Fig. 1A, while f₃ is a noisy version of f₁, derived as f₃ = f₁ + 1.5n, where n is drawn from a normal distribution $\mathcal{N}(0,1)$. Figs. 1B, 1C and 1D show in a continuous black line the PDFs derived from the data, and in a grey dash line their fitted Gaussian, in dimensions f₁, f₂ and f₃ respectively. This fitted Gaussian was derived using the Curve Fitting toolbox of Matlab™. As it can be observed, the matching of a Gaussian with an irrelevant feature is almost perfect, while those features that contain more information, like f₁ and f₂, deviate much more from a normal distribution.

Making use of these differences, an error, denoted ϕ_l, for each feature l, where $l=1,\dots ,d$, is calculated as:

(4) $${\varphi }_l={\displaystyle \frac{1}{H}\sum _{i=1}^H(p_i-g_i{)}^2,}$$

where H is the number of bins in which the range of the data is divided to estimate the PDF (p), and g is the fitted Gaussian. The number of bins in this work is set to 100 for standardization purposes. Equation (4) corresponds to the mean-squared error (MSE) between the PDF of the data over feature l and its fitted Gaussian. From these ϕ_l, the final weights b_l are calculated as:

(5) $$b_l={\displaystyle \frac{{\varphi }_l}{\sum _{l=1}^d{\varphi }_l}}$$

Therefore, combining (3) and (5), the proposed approximation, denoted ${\hat{\sigma }}^2$, is derived as:

(6) $${\hat{\sigma }}^2=\sum _{l=1}^d{b}_l{\delta }_l,$$

which gathers the distances present in the most relevant features, giving less importance to the dimensions that do not contribute to describe the structure of the data. The complete algorithm to calculate ${\hat{\sigma }}^2$ is described in Algorithm 1.

Utility metric for feature subset selection

In the manifold learning stage, a new representation E of the data based on the eigenvectors was built, which described the main structures present in the original high-dimensional data. The goal is to select a subset of the features which best approximates the data in this new representation. In the literature, this feature selection problem is formulated using a graph-based loss function and a sparse regularizer of the coefficients is used to select a subset of features, as explained in Zhu et al. (2016). The main idea of these approaches is to regress the data to its low dimensional embedding along with some sparse regularization. The use of such regularization techniques reduces overfitting and achieves dimensionality reduction. This regression is generally formulated as a least squares (LS) problem, and in many of these cases, the metric that is used for feature selection is the magnitude of their corresponding weights in the least squares solution (Cai, Zhang & He, 2010; Gui et al., 2016). However, the optimized weights do not necessarily reflect the importance of the corresponding feature as it is scaling dependent and it does not properly take interactions across features into account (Bertrand, 2018). Instead, the importance of a feature can be quantified using the increase in least-squared error (LSE) if that feature was to be removed and the weights were re-optimized. This increase in LSE, called the ‘utility’ of the feature can be efficiently computed (Bertrand, 2018) and can be used as an informative metric for a greedy backwards feature selection procedure (Bertrand, 2018; Narayanan & Bertrand, 2020; Szurley et al., 2014), as an alternative for (group-) LASSO based techniques. Under some technical conditions, a greedy selection based on this utility metric can even be shown to lead to the optimal subset (Couvreur & Bresler, 2000).

After representing the dataset using the matrix $\mathbf{E}\in {\mathbb{R}}^{N\times c}$ containing the c eigenvectors, the following LS optimization problem finds the weights p that best approximate the data X in the c-dimensional representation in E:

(7) $$J=\underset{\mathbf{P}}{min}{\displaystyle \frac{1}{N}||\mathbf{X}\mathbf{p}-\mathbf{E}|{|}_F^2}$$

where J is the cost or the LSE and ||.||_F denotes the Frobenius norm.

If X is a full rank matrix and if N > d, the LS solution $\hat{\mathbf{p}}$ of (7) is

(8) $$\hat{\mathbf{p}}={\mathbf{R}}_{\mathbf{X}\mathbf{X}}^{-1}{\mathbf{R}}_{\mathbf{X}\mathbf{E}},$$

with ${\mathbf{R}}_{\mathbf{X}\mathbf{X}}={\displaystyle \frac{1}{N}{\mathbf{X}}^T\mathbf{X}}$ and ${\mathbf{R}}_{\mathbf{X}\mathbf{E}}={\displaystyle \frac{1}{N}{\mathbf{X}}^T\mathbf{E}}$.

The goal of this feature selection method is to select the subset of s(<d) features that best represents E. This feature selection problem can be reduced to the selection of the best s(<d) columns of X which minimize (7). However, this is inherently a combinatorial problem and is computationally unfeasible to solve. Nevertheless, several greedy and approximative methods have been proposed (Gui et al., 2016; Nie, Zhu & Li, 2019; Narayanan & Bertrand, 2020). In the current work, the use of the utility metric for subset selection is proposed to select these best s columns.

The utility of a feature l of X, in an LS problem like (7), is defined as the increase in the LSE J when the column corresponding to the l-th feature in X is removed from the problem and the new optimal weight matrix, ${\hat{\mathbf{p}}}_{-l}$, is re-computed similar to (8). Consider the new LSE after the removal of feature l and the re-computation of the weight matrix ${\hat{\mathbf{p}}}_{-l}$ to be J_−l, defined as:

(9) $$J_{-l}={\displaystyle \frac{1}{N}||{\mathbf{X}}_{-l}{\hat{\mathbf{p}}}_{-l}-\mathbf{E}|{|}_F^2}$$

where X_−l denotes the matrix X with the column corresponding to l-th feature removed. Then according to the definition, the utility of feature l, U_l is:

(10) $$U_l=J_{-l}-J$$

A straightforward computation of U_l would be computationally heavy due to the fact that the computation of ${\hat{\mathbf{p}}}_l$ requires a matrix inversion of ${\mathbf{X}}_{-l}{\mathbf{X}}_{-l}^T$, which has to be repeated for each feature l.

However, it can be shown that the utility of the l-th feature of X in (10) can be computed efficiently without the explicit recomputation of ${\hat{\mathbf{p}}}_{-l}$ by using the following expression (Bertrand, 2018):

(11) $$U_l={\displaystyle \frac{1}{q_l}||{\overline{\mathbf{p}}}_l|{|}_2,}$$

where q_l is the l-th diagonal element of R⁻¹_XX and p_l is the l-th row in $\hat{\mathbf{p}}$, corresponding to the l-th feature. The mathematical proof of (11) can be found in Bertrand (2018). Note that R⁻¹_XX is already known from the computation of $\hat{\mathbf{p}}$ such that no additional matrix inversion is required.

However, since the data matrix X can contain redundant features or features that are linear combinations of each other in its columns, it cannot be guaranteed that the matrix X in (7) is full-rank. In this case, the removal of a redundant column from X will not lead to an increase in the LS cost of (7). Moreover, R⁻¹_XX, used to find the solution of (7) in (8), will not exist in this case since the matrix X is rank deficient. A similar problem appears if N < d, which can happen in case of very high-dimensional data. To overcome this problem, the definition of utility generalized to a minimum l₂-norm selection (Bertrand, 2018) is used in this work. This approach eliminates the feature yielding the smallest increase in the l₂-norm of the weight matrix when the column corresponding to that feature were to be removed and the weight matrix would be re-optimized. Moreover, minimizing the l₂-norm of the weights further reduces the risk of overfitting.

This generalization is achieved by first adding an l₂-norm penalty β to the cost function that is minimized in (7):

(12) $$J=\underset{\mathbf{p}}{min}{\displaystyle \frac{1}{2}||\mathbf{X}\mathbf{p}-\mathbf{E}|{|}_F^2+\beta ||\mathbf{p}|{|}_2^2}$$

where 0 < β μ with μ equal to the smallest non-zero eigenvalue of R_XX in order to ensure that the bias added due to the penalty term in (12) is negligible. The minimizer of (12) is:

(13) $$\hat{\mathbf{p}}={\mathbf{R}}_{\mathbf{X}\mathbf{X}\beta }^{-1}{\mathbf{R}}_{\mathbf{X}\mathbf{E}}=({\mathbf{R}}_{\mathbf{X}\mathbf{X}}+\beta \mathbf{I}{)}^{-1}{\mathbf{R}}_{\mathbf{X}\mathbf{E}}$$

It is noted that (13) reduces to R^†_XX R_XE when $\beta \to 0$, where R^†_XX denotes the Moore-Penrose pseudo-inverse. This solution corresponds to the minimum norm solution of (7) when X contains linearly dependent columns or rows. The utility U_l of the l-th column in X based on (12) is (Bertrand, 2018):

(14) $$\begin{array}{cc}{U}_l\hfill & =\left(||{\mathbf{X}}_{-l}{\hat{\mathbf{p}}}_{-l}-\mathbf{E}|{|}_2^2-||\mathbf{X}\hat{\mathbf{p}}-\mathbf{E}|{|}_2^2\right)\hfill \\ \hfill & +\beta \left(||{\hat{\mathbf{p}}}_{-l}|{|}_2^2-||\hat{\mathbf{p}}|{|}_2^2\right)\hfill \\ \hfill & =\left(J_{-l}-J\right)+\beta \left(||{\hat{\mathbf{p}}}_{-l}|{|}_2^2-||\hat{\mathbf{p}}|{|}_2^2\right)\hfill \end{array}$$

Note that if column l in X is linearly independent from the other columns, (14) closely approximates to the original utility definition in (10) as the first term dominates over the second. However, if column l is linearly dependent, the first term vanishes and the second term will dominate. In this case, the utility quantifies the increase in l₂-norm after removing the l-th feature.

To select the best s features of X, a greedy selection based on the iterative elimination of the features with the least utility is carried out. After the elimination of each feature, a re-estimation of the weights $\hat{\mathbf{p}}$ is carried out and the process of elimination is repeated, until s features remain.

Note that the value of β depends on the smallest non-zero eigenvalue of R_XX. Since R_XX has to be recomputed every time when a feature is removed, also its eigenvalues change along the way. In practice, the value of β is selected only once and fixed for the remainder of the algorithm, as smaller than the smallest non-zero eigenvalue of R_XX before any of the features are eliminated (Narayanan & Bertrand, 2020). This value of β will be smaller than all the non-zero eigenvalues of any principal submatrix of R_XX using the Cauchy’s interlace theorem (Hwang, 2004).

The summary of the utility subset selection is described in Algorithm 2. Algorithm 3 outlines the complete U2FS algorithm proposed in this paper.

As it has been stated before, one of the most remarkable aspects of the U2FS algorithm is the use of a greedy technique to solve the subset selection problem. The use of this type of method reduces the computational cost of the algorithm. This can be confirmed analyzing the computational complexity of U2FS, where the most demanding steps are the eigendecomposition of the Laplacian matrix (step 2 of Algortihm 4), which has a cost of O(N³) (Tsironis et al., 2013), and the subset selection stage in step 3 of Algorithm 4. Contrary to the state-of-the-art, the complexity of U2FS being a greedy method depends on the number of features to select. The most computationally expensive step of the subset selection in U2FS is the calculation of the matrix ${\mathbf{R}}_{\mathbf{X}\mathbf{X}}^{-1}$, which has a computational cost of O(d³). In addition, this matrix needs to be updated d − s times. This update can be done efficiently using a recursive updating equation from Bertrand (2018) with a cost of O(t²), with t the number of features remaining in the dataset, i.e., t = d − s. Since t < d, the cost for performing d − s iterations will be O((d − s)d²), which depends on the number of features s to be selected. Note that the cost of computing the least squares solution ${\hat{\mathbf{p}}}_{-l}$ for each l in (14) is eliminated using the efficient Eq. (11), bringing down the cost for computing the utility from O(t⁴) to O(t) in each iteration. This vanishes with respect to the O(d³) term (remember that t < d). Therefore, the total asymptotic complexity of U2FS is O(N³ + d³).

Simulations

A set of nonlinear toy examples typically used in clustering problems are proposed to test the different feature selection methods. In these experiments, the goal was to verify the correct selection of the original set of features. Figure 2 shows the toy examples considered¹ , which are described by features f₁ and f₂, and the final description of the datasets can be seen in Table 2.

All these problems are balanced, except for the last dataset Cres-Moon, for which the data is divided 25% to 75% between the two clusters. Five extra features in addition to the original f₁ and f₂ were added to each of the datasets in order to include redundant or irrelevant information:

• f′₁ and f′₂: random values extracted from two Pearson distributions characterized by the same higher-order statistics as f₁ and f₂ respectively.

• f′₃ and f′₄: Original f₁ and f₂ contaminated with Gaussian noise ( $\nu \mathcal{N}(0,1)$), with ν = 1.5.

• f′₅: Constant feature of value 0.

The first step in the preprocessing of the features was to standardize the data using z-score to reduce the impact of differences in scaling and noise. In order to confirm the robustness of the feature selection techniques, the methods were applied using 10-fold cross-validation on the standardized data. For each fold a training set was selected using m-medoids, setting m to 2,000 and using the centers of the clusters found as training samples. By doing so, the generalization ability of the methods can be guaranteed (Varon, Alzate & Suykens, 2015). On each of the 10 training sets, the features were selected applying the 5 methods mentioned in Table 1. For each of the methods, the number of clusters c was introduced as the number of classes presented in Table 2. Since these experiments aim to evaluate the correct selection of the features, and the original features f₁ and f₂ are known, the number of features s to be selected was set to 2.

Regarding the parameter settings within the embedding methods, the binary was obtained setting k in the kNN approach to 5. For the RBF kernel embedding, ${\sigma }_0^2$ was set to the mean of the standard deviation along each dimension, as done in Cai, Zhang & He (2010). When using ${\hat{\sigma }}^2$, its value was obtained by applying the method described in Algorithm 1.

In terms of subset selection approaches, the method based on the l₁ − norm automatically sets the value of the regularization parameter required for the LARS implementation, as described in (Cai & Chiyuan Zhang, 2020). For the utility metric, β was automatically set to the smallest non-zero eigenvalue of the matrix R_XX as described in Algorithm 2.

The performance of the algorithm is evaluated comparing the original set of features f₁ and f₂ to those selected by the algorithm. In these experiments, the evaluation of the selection results is binary: either the feature set selected is correct or not, regardless of the additional features f′_i, for i = 1,2,…,5, selected.

In Table 3 the most common results obtained in the 10 folds are shown. The utility-based approaches always obtained the same results for all 10 folds of the experiments. On the contrary, the l₁ − norm methods provided different results for different folds of the experiment. For these cases, Table 3 shows the most common feature pair for each experiment, occurring at least 3 times.

As shown in Table 3, the methods that always obtain the adequate set of features are based on utility, both with the binary weighting and with the RBF kernel and the suggested ${\hat{\sigma }}^2$. Since these results were obtained for the 10 folds, they confirm both the robustness and the consistency of the U2FS algorithm.

Benchmark datasets

Additionally, the proposed methods were evaluated using 6 well-known benchmark databases. The databases considered represent image (USPS, ORL, COIL20), audio (ISOLET) and text data (PCMAC, BASEHOCK)² , proposing examples with more samples than features, and vice versa. The description of these databases is detailed in Table 4. All these datasets are balanced, except USPS.

In these datasets, the relevant features are unknown. Therefore, the common practice in the literature to evaluate feature selectors consists of applying the algorithms, taking from 10% to 80% of the original set of features, and evaluating the accuracy of a classifier when trained and evaluated with the selected feature set (Zhu et al., 2016). The classifier used for this aim in other papers is k-Nearest Neighbors (KNN), setting the number of neighbors to 5.

These accuracy results are computed using 10-fold cross-validation to confirm the generalization capabilities of the algorithm. By setting m to 90% of the number of samples available in each benchmark dataset, m-medoids is used to select the m centroids of the clusters and use them as training set. Feature selection and the training of the KNN classifier are performed in these 9 folds of the standardized data, and the accuracy of the KNN is evaluated in the remaining 10% for testing. Exclusively for USPS, given the size of the dataset, 2,000 samples were used for training and the remaining data was used for testing. These 2,000 samples were also selected using m-medoids. Since PCMAC and BASEHOCK consist of binary data, these datasets were not standardized.

The parameters required for the binary and RBF embeddings, as well as β for the utility algorithm, are automatically set as detailed in “Discussion”.

Figure 3 shows the median accuracy obtained for each of the 5 methods. The shadows along the lines correspond to the 25 and 75 percentile of the 10 folds. As a reference, the accuracy of the classifier without using feature selection is shown in black for each of the datasets. Additionally, Fig. 4 shows the computation time for both the utility metric and the l₁ − norm applied on a binary weighting embedding. In this manner, the subset selection techniques can be evaluated regardless of the code efficiency of the embedding stage. Similarly to Fig. 3, the computation time plots show in bold the median running time for each of the subset selection techniques, and the 25 and 75 percentiles around it obtained from the 10-fold cross-validation.

The difference in the trends of the l₁ − norm and utility in terms of computation time is due to their formulation. Feature selection based on l₁ − norm regularization, solved using the LARS algorithm in this case, requires the same computation time regardless of the number of features aimed to select. All features are evaluated together, and later on, an MCFS score obtained from the regression problem is assigned to them (Cai, Zhang & He, 2010). The features with the higher scores are the ones selected. On the other hand, since the utility metric is applied in a backward greedy trend, the computation times change for different number of features selected. The lower the number of features selected compared to the original set, the higher the computation time. This is aligned with the computational complexity of the algorithm, described in “Related Work”. In spite of this, it can be seen that even the highest computation time for utility is lower than the time taken using l₁ − norm regularization. The experiments were performed with 2x Intel Xeon E5-2640 @ 2.5 GHz processors and 64GB of working memory.

Finally, the experiments in benchmark databases are extended to compare U2FS to other key algorithms in the state-of-the-art. As it was mentioned at the beginning of this section, the selected algorithms are MCFS, NDFS, and RJGSC, which represent, respectively, the precursor of U2FS, an improved version of MCFS, and an example from the class of adaptive algorithms which recursively optimize the objective function proposed. NDFS and RJGSC require the tuning of their regularization parameters, for which the indications in their corresponding articles were followed. For NDFS, the value of γ was set to 10⁸, and α and β were selected from the values {10 ^{− 6},10 ^{− 4},…,10⁶} applying grid search. The matrix F was initialized with the results of spectral clustering using all the features. For RJGSC, the results described in Zhu et al. (2016) for the BASEHOCK and PCMAC datasets are taken as a reference. In MCFS, the embedding is done using KNN and binary weighting, and the l₁ − norm is used for subset selection. U2FS, on the other hand, results from the combination of the RBF kernel with ${\hat{\sigma }}^2$ and the utility metric. Table 5 summarizes the results by showing the KNN accuracy (ACC) for 10% of the features used, and the maximum ACC achieved among the percentages of features considered, for the BASEHOCK and PCMAC datasets.