Algorithmic Stability and Generalization of an Unsupervised Feature Selection Algorithm
Year: 2021

Abstract:
Feature selection, as a vital dimension reduction technique, reduces data dimension by identifying an essential subset of input features, which can facilitate interpretable insights into learning and inference processes. Algorithmic stability is a key characteristic of an algorithm regarding its sensitivity to perturbations of input samples. In this paper, we propose an innovative unsupervised feature selection algorithm attaining this stability with provable guarantees. The architecture of our algorithm consists of a feature scorer and a feature selector. The scorer trains a neural network (NN) to globally score all the features, and the selector adopts a dependent sub-NN to locally evaluate the representation abilities for selecting features. Further, we present algorithmic stability analysis and show that our algorithm has a performance guarantee via a generalization error bound. Extensive experimental results on real-world datasets demonstrate superior generalization performance of our proposed algorithm to strong baseline methods. Also, the properties revealed by our theoretical analysis and the stability of our algorithmselected features are empirically confirmed.

heading: 1 Introduction
text: High-dimensional data is challenging due to the curse of dimensionality [7]. Dimensionality reduction is an important technique for dealing with such data, comprising two typical approaches: feature extraction and feature selection. The former, including principal component analysis (PCA) [35] and autoencoder (AE) [41, 5], is widely used in various fields such as biology [3] and computer vision [49, 18, 38]. Nonetheless, new features produced by feature extraction form a new space and, in general, do not have a direct correspondence to original features, leading to difficulty in deriving interpretable insights for the domain problems, such as biomarker identification and drug discovery. Alternatively, feature selection identifies essential features from the original feature space, providing critical interpretations and insights in many tasks [10, 45], e.g., gene functional enrichment analysis, biomarker detection, and high-throughput screening for drug discovery [22, 47, 27].
Stability is an important characteristic of an algorithm, which quantifies the sensitivity of the output to the perturbation of its training samples [12, 19]. The current stability analyses of feature selection algorithms mainly focus on similarity-based and frequency-based stability measures for practical assessment, e.g., [16, 21, 34], rather than the algorithms themselves, leaving feature selection with proper algorithmic stability [8] an unmet need. Especially, existing unsupervised feature selection algorithms use the empirical error, or its proxy, to represent generalization error but provide no theoretical guarantee of such use, leaving their ability to generalize to new data unclear.
To address these issues, in this paper we propose a novel unsupervised feature selection algorithm with a proven algorithmic stability guarantee. Our approach trains a neural network (NN) to globally score all features and a dependent sub-NN to select the features with the highest scores to reconstruct the
∗Correspondence should be addressed to: qiang.cheng@uky.edu.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
original data. Mathematically, we analyze our new algorithm and provide a generalization error upper bound, ensuring its algorithmic stability for guaranteed learning performance. Our contributions are summarized in the following:
• We propose an innovative approach for feature selection. It constructs a NN with a feature scorer to globally score all the features and a dependent sub-NN with a feature selector to locally evaluate the representation ability of highly scored features. Thus, our algorithm is capable of both globally exploring and locally excavating essential features.
• We establish performance guarantees for our algorithm, including a proven convergence rate O ( 1/ ( nmin{ √ λ1, λ1} ) + 1/n ) to ensure uniform stability and a rate
O ( 1/ (√ nmin{ √ λ1, λ1} ) + 1/ √ n )
for generalization error. Here, n is the number of training samples, and λ1 is a regularization parameter.
• We confirm the effectiveness of our proposed algorithm with extensive experiments on 10 real datasets. It achieves more competitive performance for data reconstruction and downstream classification tasks than state-of-the-art methods. Notably, the features selected by our algorithm have performance comparable to the original features. Further, the properties revealed by our theoretical analysis and the stability of our algorithm-selected features are empirically verified.
The remainder of the paper is organized as follows. We first discuss the related work, then present our proposed algorithm. Next, we show that our proposed algorithm is uniformly stable and has a guaranteed generalization bound. Finally, we conduct extensive experiments to validate our proposed new algorithm and the properties related to our proved generalization bound.

heading: 2 Related work
text: In the literature, the approaches for bounding generalization error can be generally grouped into two categories: One by controlling the complexity of hypothesis spaces and the other by focusing on the property of learning algorithms. The former usually need to define an additional complexity measure on the hypothesis space, such as the Vapnik-Chervonenkis dimension [44] and the Rademacher complexity [6]; the latter mainly include algorithmic stability [8] and robustness-based analysis [50]. In this paper, we focus on algorithmic stability to analyze the stability and generalization of our proposed new algorithm.
Algorithmic stability in learning theory. Algorithmic stability is an important tool for analyzing the generalization of algorithms in learning theory. It characterizes the sensitivity of the loss when the inputs to an algorithm are changed. Simply speaking, if an algorithm is stable, then its loss does not change significantly when the training samples are modified slightly, such as deleting and replacing a sample. Since the inception of the notion [42, 43], the concepts and properties of stability have become a significant learning theory topic. Bousquet and Elisseeff [8] have shown that stability has a direct connection with generalization in that the uniform stability of a learning algorithm implies a tight generalization error bound. Thanks to this connection, algorithmic stability has been widely used for generalization analysis of learning algorithms, including regularized least squares regression [8], multi-class support-vector-machine classification [39], multi-task learning [28], supervised autoencoders [23], and distributed learning [46]. Especially in [28], the auxiliary tasks are used for regularization, and so is the reconstruction [23]. It is noted that all these studies are on supervised learning. For unsupervised learning, stability analysis of algorithms is yet to be developed; particularly, the role of regularization in the stability bound remains elusive.
Stability studies of feature selection algorithms. A few existing studies about the stability of feature selection algorithms are mainly on stability measures for quantifying the sensitivity to input changes. In [21] a Shannon entropy-based stability measure is proposed for evaluating the stability of selected features. In [16] the stability of a feature selection algorithm is assessed by comparing the average pairwise similarity of all feature subsets obtained from different subsamples by the algorithm. The recent work [34] generalizes the requirements of stability measures into five properties and proposes a statistical estimator for stability to satisfy these properties. Although able to provide relevant assessment information about the selected features to some extent, most stability measures have no essential connection with the corresponding feature selection algorithms. In addition, they
are demonstrated empirically in general, without theoretical insights and generalization guarantees. The only exception appears to be the theoretical proof for uniform weighting stability [25]. This work builds an algorithm called feature weighting as regularized energy-based learning (FREL), proves stability, and extends to a further ensemble version. FREL is for supervised feature selection, and the uniform weighting stability does not appear to be linked with generalization. In brief, despite the wide use of algorithmic stability in analyzing many learning algorithms, its use and analysis for unsupervised feature selection algorithms are still unmet needs.
Unsupervised feature selection algorithms. Many unsupervised feature selection algorithms have been proposed in recent years to meet the demand for effectively handling large-scale and highdimensional data. They can be generally classified into four categories: filter, wrapper, embedder, and hybrid approaches [2, 24]. In this paper, we will mainly compare our algorithm with ten typical unsupervised feature selection algorithms, including Laplacian score (LS) [17], principal feature analysis (PFA) [30], SPEC [52], multi-cluster feature selection (MCFS) [9], unsupervised discriminative feature selection (UDFS) [51], nonnegative discriminative feature selection (NDFS) [26], Autoencoder feature selector (AEFS) [15], agnostic feature selection with slack variables (AgnoS-S) [11],2 graph-based infinite feature selection (Inf-FS) [40], and concrete autoencoders (CAE) [1]. LS, SPEC, and Inf-FS are filter approaches; PFA can be classified as a wrapper approach; the other six are embedded algorithms. Among them, AEFS, AgnoS-S, and CAE are AE-based feature selection methods: AEFS combines AE regression and `2,1 regularization on the weights of the encoder to obtain a subset of useful features; AgnoS-S adopts AE with the `1 norm on slack variables in the first layer of AE to implement feature selection; CAE replaces the first hidden layer of AE with a concrete selector layer [31], and then it selects the features with a high probability of connection to the nodes of the concrete selection layer. While frequently used, these methods have no theoretical analysis of their algorithmic stability and leave their generalization abilities to new data unclear. Indeed, they generally lack such needed abilities, as demonstrated in our experiments (see Table 3).

heading: 3 Notations and preliminary
text: Let n, m, k, and d be the numbers of training samples, features, selected features, and reduced dimensions. Let X ∈ Rn×m be a sample matrix, and XT be its transposition. A lowercase capital letter like wm denotes a vector; Diag(wm) represents a diagonal matrix with the diagonal wm. wmaxkm stands for the operation to keep the k largest entries of wm while making the other entries 0. Let φ : Rm → Rm be an element-wise operation. ‖ · ‖F, ‖ · ‖2, ‖ · ‖1, and ‖ · ‖∞ respectively denote the Frobenius, `2, `1, and infinity norms. The sample space of data is denoted as X . We assume all samples are bounded, i.e., ∀x ∈ X , ∃κ1 > 0, such that ‖x‖2 6 κ1. Besides, ∀x, x′ ∈ [0,∞), for a function f , ∆t(f(x), x′) , f(x)− f(x+ t∆x), where ∆x = x′ − x and t ∈ [0, 1].
Let S , {xi ∈ X , i = 1, 2, . . . , n} be a finite set of training samples which are independently and identically distributed according to an unknown distribution P. We denote the sets after removing and replacing the i-th element from S by S\i , {x1, . . . , xi−1, xi+1, . . . , xn} and Si , {x1, . . . , xi−1, x′i, xi+1, . . . , xn}, respectively. Let A be an algorithm, and AS be a function picked by A from the hypothesis space H based on S. A is assumed to be deterministic and symmetric with respect to S. Let ` : X × X → [0,+∞) be a loss function. The generalization error is
L(A, S) , Ex [`(AS , x)] = ∫ X `(AS , x)dP,
and the empirical error is
Lemp(A, S) , 1
n n∑ i=1 `(AS , xi),
where xi ∈ S, i = 1, 2, . . . , n. To study the stability of A, the leave-one-out error is used,
Lloo(A, S) , 1
n n∑ i=1 `(AS\i , xi).
2In [11], three variants of AE-based agnostic feature selection algorithms are proposed. AgnoS-S is generally the best among the three; thus, in this paper we will compare our algorithm with AgnoS-S.
By the empirical and leave-one-out errors, uniform stability is defined as follows: Definition 1 (Uniform Stability [8]). An algorithm A has uniform stability β with respect to ` if ∀S ∈ Xn, i = 1, 2, . . . , n, ‖`(AS , ·)− `(AS\i , ·)‖∞ 6 β, where β is a function of n. Generally, for a uniformly stable algorithm, β decreases as O(1/n).

heading: 4 Method
text: Now we will present our novel AE-based feature selection algorithm.
Revisit of AE. For AE, we formalize it as follows:
min f,g ‖X− f(g(X))‖2F,
where g is an encoder, and f is a decoder. g(X) embeds the input data into a latent space Rn×d, where d denotes the dimension of the bottleneck layer.
Formalization of unsupervised feature selection. The goal of feature selection is to identify a subset of important features in the original feature space, and it can be formalized as follows:
min S(k),H
‖H(XS(k))−X‖2F, (1)
where S(k) denotes a subset of k features of X with k < m, XS(k) is the resulting dataset by restricting X to the selected subset. We use H to represent a mapping on the k-dimensional space, and it selects a subset of the features of X to preserve the information of X as much as possible. The optimization problem of evaluating the subset of features is typically NP-hard [32, 14]. This paper will develop an efficient algorithm to effectively approximate the solution of (1).
AE-based new algorithm for feature selection. Our model is constructed as follows:
min WI,f,g
‖X− f(g(X(Φ(WI)maxk)))‖2F + λ1‖X− f(g(X(Φ(WI))))‖2F, (2)
where Φ(WI) , Diag(φ(wm)) ∈ Rm×m, Φ(WI)maxk , Diag((φ(wm))maxk) ∈ Rm×m, and λ1 is a regularization parameter. We use NN for the optimization of (2), and the network architecture is shown in Supplementary Figure 1. For the convenience of later discussions, ∀x ∈ X , let `selec(Φ(WI)maxk , x) , ‖x − f(g(x(Φ(WI)maxk)))‖22, and `score(Φ(WI), x) , ‖x− f(g(x(Φ(WI))))‖22. The operators Φ(WI) and Φ(WI)maxk are called feature scorer and feature selector, respectively.3 Φ(WI) and Φ(WI)maxk will iteratively interact through the NN and sub-NN. Specifically, (2) can be regarded as an extension of a common NN, with the main difference being the sub-NN, i.e., the first term of (2). During training, Φ(WI)maxk will select the top-k features (from the scorer weights, i.e., Φ(WI)) by ensuring these selected features to well reconstruct the original input X with the NN. On the other hand, the NN, as a constraint, will ensure the selected features for Φ(WI)maxk to have the largest importances globally. After training, we obtain the trained Φ(WI)maxk , and then we use it to make feature selection on new samples during testing. Taking MNIST-Fashion, COIL-20, and USPS as examples, for k = 50 and Φmaxk = (W2I )
maxk , we visualize the feature selection and reconstruction results on MNIST demonstrated in Figure 1, and we give the selected features on original samples of USPS in Figure 2.

heading: 5 Algorithmic stability and generalization
text: As previously reviewed, existing feature selection methods have not considered algorithmic stability. Generally, they lack this property or generalization ability, as demonstrated in our experiments. Moreover, there has been no clue about the role of regulation in uniform stability. Our algorithm is built to address these issues and, here, it will be shown to possess these desirable properties. For (2),
3For simplicity, we will use the shorthand notations Φ and Φmaxk to denote Φ(WI) and Φ(WI)maxk .
the first term implies that the original samples are well represented by the selected k features. The second term scores all the features by their reconstruction ability and ensures that the selected features have high scores among all features. It will be used as a regularizer to promote the stability of the feature selection. It requires that all features be explored globally for high representation ability before being further excavated locally by the selector. The role of the scorer will be empirically demonstrated in the experiments (see Figure 3 (b)).
It is worth noting that, in contrast to existing models, such as NNs, which require nonlinear activation functions to capture the nonlinearity, the first term in (2) naturally facilitates a nonlinear model even when the scorer is a linear model. Such a nonlinearity intrinsically comes from ranking the feature weights and selecting the k informative features during the iterative optimization. For simplicity, in the following discussion, we only consider the linear version of f and g in (2), that is, g(X) = XWE, WE ∈ Rm×d, and f(g(X)) = (g(X))WD, WD ∈ Rd×m. The experimental section will show that such a concise linear version already achieves superior performance.
Formalization of stability and generalization. Let (Φmaxk∗ ,Φ∗) and ((Φ \j ∗ ) maxk ,Φ \j ∗ ) correspond to the optimal feature selectors and feature scorers for the errors
RLemp(Φ, S) , 1
n n∑ i=1 `selec(Φmaxk , xi) + λ1` score(Φ, xi)
and
RL\jemp(Φ, S) , 1
n n∑ i=1,i6=j `selec(Φmaxk , xi) + λ1` score(Φ, xi),
where j ∈ {1, . . . , n}. Our goal is to bound `selec(Φmaxk∗ , ·)−`selec((Φ \j ∗ )
maxk , ·), which will be used for the analysis of uniform stability. Further, we intend to boundLselec(Φmaxk∗ , S)−Lselecemp (Φmaxk∗ , S), or Lselec(Φmaxk∗ , S)− Lselecloo (Φmaxk∗ , S), to be used in the analysis of generalization errors.
Assumptions. To determine the stability and generalization bounds, the following assumptions will be needed: Assumption 1. ∃κ2 > 0, such that, ∀x ∈ X , ‖(x(Φmaxk∗ − (Φ \j ∗ )
maxk)WE)WD‖2 6 κ2‖(x(Φ∗ − Φ \j ∗ )WE)WD‖2.
This assumption is reasonable, because the null space of Φ∗−Φ\j∗ is smaller than Φmaxk∗ −(Φ \j ∗ )
maxk ; that is, we generally have
r (((( Φmaxk∗ − (Φ \j ∗ ) maxk ) WE ) WD )T) 6 r ((( Φ∗ − Φ\j∗ ) WEWD )T) ,
where r(·) denotes the rank of a matrix. So, it is more likely that the column vector ((Φ∗ − Φ \j ∗ )WEWD)
TxT will have m non-zero elements. Consequently, it is more likely that the row vector (x(Φ∗ − Φ\j∗ )WE)WD will have m non-zero elements.
Assumption 2. ∃Z = {z1, z2, . . . , zu} ⊂ S, such that, ∀x ∈ X , x can be reconstructed by Z, i.e., x = ∑u i=1 αizi + η, where η is a small reconstruction error that satisfies ‖η‖2 6 κ3/n, and αi is a
scalar that satisfies √∑u
i=1 α 2 i 6 κ4, with κ3 and κ4 being positive constants.
This assumption is mild and implies that any sample x can be reconstructed by Z with a small construction error that decreases with n.4 Theoretically speaking, when the samples all reside on a manifold, it obviously holds. Practically, it is related to the self expressiveness of samples exploited in many subspace learning algorithms, e.g., [29, 36, 37].
Assumption 3. Let Lscore,\jemp (Φ, S) , (1/n) ∑n i=1,i6=j `
score(Φ, xi) and Lscoreemp,z(Φ) , (1/u) ∑u i=1 ` score(Φ, zi). ∀t ∈ [0, 1],
∆t(Lscoreemp,z(Φ∗),Φ \j ∗ ) + ∆ t(Lscoreemp,z(Φ \j ∗ ),Φ∗)
6 ∆t(L
score,\j emp (Φ∗, S),Φ \j ∗ ) + ∆ t(L score,\j emp (Φ \j ∗ , S),Φ∗)
(n− 1)/n .
(3)
The implication of this assumption is clear: it requires that, at Φ∗ and Φ \j ∗ , the perturbations of Lscoreemp,z can be controlled by those of Lscoreemp , and the left-hand side of (3) is related to Z in Assumption 2. [23] makes a similar assumption; however, we enhance the factor t into n/(n − 1) on the right-hand side of (3), which makes Assumption 3 much weaker and more reasonable than that in [23], since Lscoreemp,z(Φ) = (n/(n−1))L score,\j emp (Φ, S) whenZ = S\j . Further, the perturbations of Lscoreemp (Φ∗, S) and Lscore,\jemp (Φ∗) are implicitly related to λ1 (see Proposition 1 below).
Proposition 1. Let ∆Φmaxk , (Φ\j∗ )maxk − Φmaxk∗ . ∀t ∈ [0, 1], the following inequality holds:
∆t(Lscoreemp (Φ∗, S),Φ \j ∗ ) + ∆ t(Lscore,\jemp (Φ∗, S),Φ \j ∗ )
6
( tκ21 (∥∥∥(2(Φ\j∗ )maxk − t∆Φmaxk)WEWD∥∥∥ 2 + 2 ) ‖((∆Φmaxk) WE) WD‖2 ) λ1 .
Proposition 1 reveals the relationship of the regularization represented by λ1 with the perturbation of Lscoreemp (Φ∗, S) and L score,\j emp (Φ∗, S).
As discussed above, we will mainly verify Assumption 2. For this purpose, we develop a coresubspace learning procedure. More details about the verification are given in Supplementary Material.
Uniform stability. We establish a bound on uniform stability of our algorithm for (2). Theorem 1 (Uniform Stability). Under Assumptions 1, 2, and 3, we have, ∀n > 2,∥∥∥`selec (Φmaxk∗ , ·)− `selec ((Φ\j∗ )maxk , ·)∥∥∥∞ = O ( 1 nmin{ √ λ1, λ1} + 1 n ) . (4)
This theorem quantifies the insensitivity of our algorithm to the perturbation of its input, and the optimal feature selector Φmaxk∗ would not change significantly with the change of one sample. The upper bound in (4) is mainly about the approximate behavior of n, and it is theoretically more meaningful when λ1 is much larger than 1/n. If λ1 is very large, i.e., with over-regularization, then essentially the second term of (2) works, which will make the feature selector underfit training samples and incur large test error. Theoretically, it may still be sufficient for generalization, and the bound in (4) will reduce to O (1/n); in practice, however, the corresponding validation error should be large. Thus, such a large λ1 would never be practically chosen. How to analytically determine the optimal value of λ1 is out of the scope of this paper and will be an interesting future research topic. In this paper, we empirically choose a positive value of λ1 by using cross-validation (see Experiments for details). If λ1 is set to 0, then the second term of (2) will vanish; in this case, the stability with only the first term of (2) warrants further study in the future.
4Liu et al [28] and Le et al [23] have made a similar assumption.
Generalization error bound. Now we will bound Lselecemp (Φmaxk∗ , S)− Lselec (Φmaxk∗ , S). Theorem 2 (Generalization Error). ∃κ5 > 0 and δ ∈ (0, 1), ∀x ∈ X and S, as long as `(AS , x) 6 κ5, the following inequality holds with probability at least 1− δ,
Lselec (Φmaxk∗ , S)− Lselecemp (Φmaxk∗ , S) = O
 √ ln ( 1 δ ) √ nmin{ √ λ1, λ1} + √ ln ( 1 δ ) n  . This theorem shows that, besides the uniform stability bound in Theorem 1, there is an upper bound for the generalization error of our proposed algorithm. By this theorem, as long as our feature selection algorithm is stable, the empirical error has a proven guarantee to approximate generalization error. The convergence rate of generalization error is O ( 1/ (√ nmin{ √ λ1, λ1} ) + 1/ √ n ) . It is noted that, when λ1 > 1, the convergence rate is O (1/ √ n), which is different from the rate under
the `2 regularization; when λ1 6 1, the convergence rate is O (1/ ( √ nλ1)), the same as that under the `2 regularization [8]. When the reconstruction error from the feature scorer is made small with a large λ1 > 1, the scorer will play a more important role than the selector, and the selector might be underfitted (see Figure 3 (b)); and vice versa for λ1 6 1. Additionally, if we fix λ1, the bound will reduce to O (1/ √ n), which decays similarly to the bound under the `2 regularization [8] when λ1 is fixed.
The proofs of Proposition 1, Theorems 1 and 2, the bound of Lselec (Φmaxk∗ , S)− Lselecloo (Φmaxk∗ , S), and more discussions are provided in Supplementary Material.

heading: 6 Experiments
text: In this section, we will perform extensive experiments to validate our new algorithm.
Datasets to be used. The benchmarking datasets and their statistics are summarized in Table 1.5 Following CAE [1] and considering the long runtime of UDFS, for dataset 5, we randomly choose 6, 000 samples from the training set for training and validating and 4, 000 samples from the test set for testing. We then randomly split 6, 000 samples into training and validation sets by a ratio of 90:10. For other datasets, we randomly split the samples into training, validation, and test sets by a ratio of 72:8:20, and we tune hyperparameters on the validation set. More details about these datasets are provided in Supplementary Material.
Design of experiments. In all experiments, for our model, we use the same linear decoder as that in [1]. For the encoder, for a fair comparison, our algorithm uses the structure from an AE having one hidden layer with the same number of neurons as other AE-based baselines in comparison. We use the linear activation function for the encoder in the same way as CAE and AEFS in [1]. For feature selection, the first layer of the encoder in our NN uses the slack variables in the same way as AgnoS-S in [11]. We set the maximum number of epochs to 200. We initialize the weights of the feature selection layer by sampling uniformly from U[0.999999, 0.9999999] and the other layers with the Xavier normal initializer.6 We adopt the Adam optimizer [20] with a learning rate of 0.001. We set λ1 to 1/27.7 We take k = 10 for dataset 1 and k = 50 for datasets 2-6 following CAE [1], and
5Datasets 1 and 4 are downloaded from http://archive.ics.uci.edu/ml/datasets/. Datasets 6-10 are from the scikit-feature feature selection repository [24].
6This distribution makes initial values of the weights in the selection layer close to 1 but still different to break the potential ties.
7We tune λ1 by searching in { 2, 1/20, . . . , 1/210 }
on the validation set of MNIST-Fashion, then choose the optimal one and use it on other datasets.
k = 64 for high-dimensional datasets 7-10. The dimension of the latent space is consistently set to k. The main codes related to our proposed algorithm are publicly available,8 and the implementation details of baseline algorithms are provided in Supplementary Material.
Two metrics are used for evaluating the performance of algorithms: 1) reconstruction error, which is measured in mean squared error (MSE); 2) accuracy, which is measured by passing the selected features to a downstream classifier as a viable means to benchmark the quality of selected features. For a fair comparison, following CAE [1], after selecting the features, we train an ordinary linear regression model to reconstruct the original features, and the resulting linear reconstruction error is used as metric 1);9 for metric 2) we adopt extremely randomized trees [13] as the classifier.
Results on 10 datasets. Our experimental results on reconstruction and classification with the selected features are displayed in Tables 2 and 3,10 and |WI| and W2I are shorthands for feature selectors Φmaxk = |WI|maxk and Φmaxk = (W2I )maxk . We followed the way of CAE [1] to split samples and to report the final results on the hold-out test set. For a fair comparison, we directly adopt the published results for reconstruction and classification by LS, AEFS, UDFS, MCFS, PFA, and CAE on datasets 1-5 from [1]. From Table 2, it is seen that our algorithm gives smaller linear reconstruction errors than baseline methods on majority datasets, indicating a stronger ability to select a subset of representative features. From Table 3, it is evident that our algorithm exhibits almost consistently superior performance in the downstream classification task on diverse datasets.
Further, we compare the behavior of our algorithm with contemporary algorithms over different k. By varying k on ISOLET we obtain the corresponding linear reconstruction errors and classification accuracy rates as the outputs from different algorithms. We plot the linear reconstruction errors in MSE and classification accuracy rates in Supplementary Figure 8. It can be observed that our algorithm demonstrates almost consistently better and more stable performance than feature selection algorithms in comparison.
8They can be found at https://github.com/xinxingwu-uk/UFS 9Here, the reconstruction error denotes the error from the first term of (2).
10The “\” mark denotes the case with prohibitive running time, where the algorithm ran for more than a week without getting a result and thus was stopped.
Also, we take Φmaxk = |WI|maxk , and we compare the variation of classification accuracy with the reduction of original features. The results are shown in Supplementary Table 1 and Supplementary
Figure 9. It is observed that all the reduction in the number of features is more than 80%, while the corresponding reduction in classification accuracy is less than 10% with the exception of ISOLET (which is 11.2%). It shows that our algorithm can effectively reduce the number of features but still maintain the classification performance comparable to the original data.

heading: 7 Discussion
text: In this section, we will empirically verify the properties related to uniform stability bound and generalization bound in Theorems 1 and 2. Further, we will empirically discuss the algorithmic stability and the stability of selected features and also analyze the time complexity of our algorithm (2).
Effect of n. We plot the curves of error difference Lselec (Φmaxk∗ , S)− Lselecemp (Φmaxk∗ , S) and test error versus n in Figure 3 (a). It is observed that n has a direct effect on the error difference and test error: 1) The more training samples, generally the smaller the test error. This observation verifies that small training sets cause overfitting and, by increasing the training set size, the generalization of our algorithm becomes improved with the test error decreased. 2) The error difference decreases in n. This phenomenon is well in line with the theoretical result in Theorem 2.
Effect of λ1. We plot the curves of the error difference and test error versus λ1 in Figure 3 (b). Two observations can be made: 1) The larger λ1, the smaller the error difference. Such an experimental result aligns well with Theorem 2. 2) The test error decreases at first and then increases with λ1. The reason may be as follows: by increasing λ1 from zero to about 0.01, the feature scorer will play an increasingly important role and facilitate finding appropriate scores for features to reduce the reconstruction error, which helps decrease the test error. The regularization role of the scorer is evident: when λ1 increases from 0 to a small value, the test error decreases. But if λ1 increases too much, it would over-emphasize the scorer and the minimization of the reconstruction error, making the feature selector underfit training samples and thus incurring larger test errors. These empirical curves offer informative clues for understanding the role of regularization in the generalization bound and choosing a proper λ1 in our experiments.
Effect of k. We plot the curves of the error difference and test error versus k in Figure 3 (c). It is seen that the larger k, the larger the error difference. This observation is in line with Theorem 2 (see Supplementary Material for the detailed expression of upper bound). As discussed below this theorem, the generalization error bound of our algorithm has a similar convergence rate to that under the `2 regularization. To see how the selection affects the error difference, we plot the curves of Frobenius norms of WEWD and W0WEWD in Figure 3 (d).11 Note that they have a similar increase tendency to the error difference. The reason may be that, by selecting more features, Lselec (Φmaxk∗ , S) and Lselecemp (Φ maxk ∗ , S) would depend on more variables and be more easily affected by the perturbations from the variables; as a result, the feature selector would be more difficult to fit. From Figure 3 (c), we can also observe that the larger k, the smaller the test error, which is sensible because more selected features mean less loss of original data and the reconstruction would be better. In addition, from Figure 3 (a)-(c), different Φ∗ scorers, such as |WI| and W2I , have similar behaviors, which reflects the tendency implied in Theorem 2.
Algorithmic stability analysis. Adopting the same experiment design on MNIST-Fashion in Section 6, we vary n from 3, 000 to 6, 000 with a step size of 1, 000 to obtain different S; meanwhile, we delete a sample for each S to get the corresponding S\i and then calculate the left-hand side of (4) on the testing set for these trained models. From the plots in Figure 3 (e)-(f), it is seen that with the increase of n, the curve of uniform stability bound presents a downward tendency, which is consistent with our theoretical analysis.
Overall, the above four experiments for interpreting the generalization bound are related to the training sample size n, the regularization parameter λ1, and the number of selected features k, and the stability bound, showing different aspects or behaviors of our algorithm. Taken together, these experiments verify the properties revealed by Theorems 1 and 2.
11Here, W0 denotes the entries of W are zeros except k ones on its diagonal.
Stability analysis of selected features. We empirically analyze the stability of features selected by (2). We randomly split the samples of MNIST-Fashion into the training and testing sets and then use (2) to perform feature selection. We repeat this procedure 10 times with different random seeds and plot the selection results in Figure 3 (g)-(h). Note that the selected features almost overlap for different splits and are stable. More results are provided in Supplementary Material.
Computational complexity. Experimentally, the computational time of our algorithm (2) is about twice that of AE. Our algorithm has only an additional sub-NN compared to AE and shares parameters with the NN. Also, the fitting error term of sub-NN is quadratic and similar to that of AE. Therefore, the overall computational complexity of (2) is of the same order as AE.
Ethical statement. This paper focuses on feature selection, which is a dimensionality reduction approach. To our best knowledge, there are no ethical issues and negative societal impacts of the proposed technique in this paper.

heading: 8 Conclusions
text: In this paper, we propose an innovative unsupervised feature selection algorithm with provable performance guarantees, which consists of a feature scorer and a feature selector. Theoretically, we prove uniform stability and provide the generalization error upper bound for our algorithm. Empirically, we show that our algorithm achieves performance better than the contemporary algorithms on various real-world datasets; additionally, the selected features by the proposed algorithm show comparable performance to the original features. Moreover, we experimentally verify the properties revealed by our theoretical analysis with respect to sample size, regularization levels, the number of selected features, and uniform stability. Additionally, this new algorithm may be applied to other tasks, such as selecting important patterns of image data beyond pixels, which we will extend in future work.

heading: 9 Funding
text: This work was partially supported by the NIH grants R21AG070909, R56NS117587, R01HD101508, and ARO W911NF-17-1-0040.

heading: Acknowledgment
text: We sincerely thank the anonymous reviewers and AC for their valuable comments.

Review summaries:
This paper proposes an unsupervised feature selection method. The paper provides theoretical analysis of selection performance when the underlying dataset is perturbed. It provides experimental results on multiple datasets, comparing the performance of the proposed approach with that of the methods in the literature.

The authors address the problem of unsupervised feature selection. They combine an autoencoder reconstruction loss with an additional term which only uses the k most informative features to reconstruct the input data. The k most informative features are identified by sorting the absolute values of a diagonal matrix whose weights are multiplied by the input features. The authors prove that the method is algorithmically stable and find its generalization error bound. The method demonstrates superior capabilities compared with several feature selection baselines.

The authors propose a novel methodology for unsupervised feature selection. Their procedure is based on simultaneously training a feature importance scorer neural network and a feature selector neural network; both are embedded inside an encoder in an autoencoder. The overall network is trained by optimizing a regularized estimator, where the regularization parameter induces some algorithmic stability. The authors provide theoretical guarantees of the estimator in the form of upper-bounds for the algorithmic stability and generalization. Further, they demonstrate the effectiveness of their approach on multiple real datasets and empirically understand various features of their algorithm (e.g. sensitivity to choices of parameters) via demonstrations on synthetic data.