Topological Feature Selection

We prove TFS’ effectiveness by running a battery of tests on $16$ benchmark datasets (all in a tabular format) belonging to $7$ different application domains (i.e. text data, face images data, hand written images data, biological data, digits data, spoken letters data and artificial data) (see Appendix A for further details). The data pre-processing pipeline consists of $3$ different steps: (i) data reading and format unification; (ii) training/validation/test splitting; (iii) constant features pruning. The first step allows to read data coming from different sources and unify the format. The second step organises each dataset into a training, validation and a test set. In the third step, non-informative, constant covariates are detected on the training set and permanently removed from the training, validation and the test set. In Appendix B, we report datasets’ specifics after the pre-processing step. We remark that most of the considered datasets are not affected by the constant features filtering step. Datasets are roughly balanced and the raw labels’ distribution is preserved through the execution of stratification (Zeng & Martinez, 2000) during the training/validation/test splitting stage.

Information Filtering Networks (IFNs) (Mantegna, 1999; Aste et al., 2005; Tumminello et al., 2005; Barfuss et al., 2016; Massara et al., 2017) are an effective tool to represent and model dependency structures among variables characterising complex systems. In the past two decades, they have been extensively used in finance (Procacci & Aste, 2022; Briola et al., 2022; Wang & Aste, 2022a; Seabrook et al., 2022; Wang & Aste, 2022b; Vidal-Tomás et al., 2023), psychology (Christensen et al., 2018; Christensen, 2018), medicine (Hutter et al., 2019; Danoff et al., 2021) and biology (Song et al., 2008, 2012). Sometimes IFNs are also referred to as Correlation Networks (CNs). Such an association is, however, inaccurate. Indeed, the two methodologies slightly differ, with CNs being normally obtained by imposing a threshold that retains only the largest correlations among variables of the system, while IFNs being constructed imposing additional topological constraints (e.g. being a tree or a planar graph) and optimising specific global properties (e.g. likelihood) (Aste, 2022). Both methodologies end in the determination of a sparse adjacency matrix, A, representing relations among variables in the system with the fundamental difference that the former approach generates a disconnected graph, while the latter guarantees the connectedness. Based on the nature of the relationships to be modelled (i.e. linear, non-linear), one can choose different metrics to build the adjacency matrix A. In most cases, A is built on an arbitrary similarity matrix $\hat{\textbf{C}}$, which often corresponds to a correlation matrix. From a network science perspective, $\hat{\textbf{C}}$ can be considered as a fully connected graph where each variable of the system is represented as a node, and each pair of variables is joined by a weighted and undirected edge representing their similarity. Historically, the main IFNs were the Minimum Spanning Tree (MST) (Papadimitriou & Steiglitz, 1998; Mantegna, 1999) and the Planar Maximally Filtered Graph (PMFG) (Tumminello et al., 2005). MSTs are a class of networks connecting all the vertices without forming cycles (i.e. closed paths of at least three nodes) while retaining the network’s representation as simple as possible (i.e. representing only relevant relations among variables characterising the system under analysis) (Briola & Aste, 2022). Prim’s algorithm for MST construction sorts all edges’ weights (i.e. similarities) in descending order and adds the largest possible edge weight among two nodes in an iterative way. The resulting network has $n-1$ edges and retains only the most significant connections, assuring, at the same time, that connectedness’ property is fulfilled (Christensen et al., 2018). Despite being a powerful method to capture meaningful relationships in network structures describing complex systems, MST presents some aspects that can be unsatisfactory. Paradoxically, the main limit is represented by its tree structure (i.e. it cannot contain cycles) which does not allow to represent direct relationships among more than two variables showing strong similarity. The introduction of the Planar Maximally Filtered Graph (PMFG) (Tumminello et al., 2005) overcomes such a shortcoming. Similarly to MST, also the PMFG algorithm sorts edge weights in descending order, incrementally adding the largest ones while imposing planarity (Christensen et al., 2018). A graph is planar if it can be embedded in a sphere without edges crossing. Thanks to this, the same powerful filtering properties of the MST are maintained, and, at the same time, extra links, cycles and cliques (i.e. complete subgraphs) are added in a controlled manner. The resulting network has $3n-6$ edges and is composed of three- and four-nodes cliques. A nested hierarchy emerges from these cliques (Song et al., 2011): dimensionality is reduced in a deterministic manner while local information and the global hierarchical structure of the original network are retained. The PMFG represents a substantial step forward compared to the MST. However, it still presents two limits: (i) it is computationally costly and (ii) it is a non-chordal graph. A graph is said to be chordal if all cycles made of four or more vertices have a chord, reducing the cycle to a set of triangles. A chord is defined as an edge that is not part of the cycle but connects two vertices of the cycle itself. The advantage of chordal graphs is that they fulfill the independence assumptions of Markov (i.e., bidirectional or undirected relations) and Bayesian (i.e., directional relations) networks (Koller & Friedman, 2009; Christensen et al., 2018). The Triangulated Maximally Filtered Graph (TMFG) (Massara et al., 2017) has been explicitly designed to be a chordal graph while retaining the strengths of PMFG. The building process of TMFG (see Appendix C for further details) is based on a simple topological move that preserves planarity: it adds one node to the centre of three-nodes cliques by using a score function that maximises the sum of the weights of the three edges connecting the existing vertices. This addition transforms three-nodes cliques (i.e. triangles) into four-nodes cliques (i.e. tetrahedrons) characterised by a chord that forms two triangles and generates a chordal network (Christensen et al., 2018). Also in this case, the resulting network has $3n-6$ edges and is composed of three- and four-nodes cliques. TMFG has two main advantages compared to PMFG: (i) it can be used to generate sparse probabilistic models as a form of topological regularization (Aste, 2022) and (ii) it is computationally more efficient. On the other hand, the two main limitations of chordal networks are that (i) they may add unnecessary edges to satisfy the property of chordality and (ii) their building cost can vary based on the chosen optimization function.

Topological Feature Selection (TFS) algorithm is a graph-based filter method to perform feature selection in an unsupervised manner. Given a set of features $F=\{f_{1},\dots,f_{n}\}$, where $n\gg 1$ is the cardinality, we build the adjacency matrix A of the corresponding TMFG based on one of the following three metrics: (i) the Pearson’s estimator of the correlation coefficient, (ii) the Spearman’s rank correlation coefficient and (iii) the Energy coefficient (i.e. the weighted combination of two pairwise measures described later in this Section). Depending on the metric’s formulation, it is possible to capture different kinds of interactions among covariates (e.g. linear or non-linear interactions).

The Pearson’s estimator of the correlation coefficient for the two covariates $f_{i}$ and $f_{j}$, is defined as:

where $S$ is the sample size, $f_{i,s}$ and $f_{j,s}$ are two sample points indexed with $s$, $\hat{\mu}$ is the sample mean and $\hat{\sigma}$ is the sample standard deviation. By definition, $r_{f_{i},f_{j}}$ has values between $-1$ (meaning that the two features are completely, linearly anti-correlated), and $+1$ (meaning that the two features are completely, linearly correlated). When $r_{f_{i},f_{j}}=0$, the two covariates are said to be uncorrelated. The Person’s estimator of the correlation coefficient heavily depends on the distribution of the underlying data and may be influenced by outliers. In addition to this, it only captures linear dependency among variables, restricting its applicability to real-world problems where non-linear interactions are often relevant (Shirokikh et al., 2013).

The Spearman’s rank correlation coefficient is based on the concept of ‘variables ranking’. Ranking a variable means mapping its realizations to an integer number that describes their positions in an ordered set. Considering a variable with cardinality $|s|$, this means assigning 1 to the realization with the highest value and $|s|$ to the realization with the lowest value.

The Spearman’s rank correlation coefficient for the two covariates $f_{i}$ and $f_{j}$, is defined as the Pearson’s correlation between the ranks of the variables:

where $S$ is the sample size, $R_{f_{i,s}}$ and $R_{f_{j,s}}$ are the ranks of the two sample points indexed with $s$, $\hat{\mu}$ is the sample mean for $R_{f_{i}}$ and $R_{f_{j}}$ and $\hat{\sigma}$ is the sample standard deviation for $R_{f_{i}}$ and $R_{f_{j}}$. If there are no repeated data samples, a perfect Spearman’s rank correlation (i.e. $r_{s_{f_{i},f_{j}}}=1$ or $r_{s_{f_{i},f_{j}}}=-1$ occurs when each of the features is a perfect monotone function of the other). Spearman’s rank correlation technique is distribution-free and allows to capture monotonic, but not necessarily linear, relationships among variables (Shirokikh et al., 2013).

The Energy coefficient is a metric introduced by (Roffo et al., 2015), and it is used as the primary benchmark for comparison between our method and the current state-of-the-art. It is a weighted combination of two different pairwise measures defined as follows:

where $E_{f_{i},f_{j}}=\max(\hat{\sigma}_{f_{i}},\hat{\sigma}_{f_{j}})$, with $\hat{\sigma}$ being the sample standard deviation computed on features $f_{i}$ and $f_{j}$ normalized to the range $[0,1]$, $\rho_{f_{i},f_{j}}=1-|r_{s_{f_{i},f_{j}}}|$ and $\alpha$ is a threshold value with a value $\in[0,1]$. $\phi_{f_{i},f_{j}}\in[0,1]$ analyses two features distributions (i.e. $f_{i}$ and $f_{j}$) taking into account both their maximal dispersion (i.e. standard deviation) and their level of uncorrelation. Computing $\phi_{f_{i},f_{j}}$ for all the features in $F$ in a pairwise manner, one can define a matrix which is $n\times n$ symmetric and completely characterized by $n(n-1)/2$ coefficients. For simplicity, we refer to this matrix as $\hat{\textbf{C}}$ too.

Once one of the above mentioned metrics is chosen and $\hat{\textbf{C}}$ is computed, TFS applies the standard TMFG algorithm defined in Appendix C on the corresponding fully connected graph, creating a sparse chordal network which is able (i) to retain useful relationships among features, (ii) prune the weakest ones, and (iii) express a significant level of information flow among input variables.

The last step toward the selection of the most relevant features, is represented by the choice of the right nodes inside the TMFG. In this sense, multiple approaches of increasing complexity can be formulated. In this paper, which is a foundational one, we study the relative position of the nodes in the network computing their degree centrality. Degree centrality is the simplest and least computationally intensive measure of centrality. Typically, all the other centrality measures are strictly related (Lee, 2006; Valente et al., 2008). Given the sparse adjacency matrix A representing the TMFG, degree centrality of a node $v$ is denoted as $\deg(v)$ and represents the number of neighbours (i.e. how many edges a node has) of $v$ as follows:

where $n$ is the cardinality of $F$ and $f_{v}$ and $f_{w}$ are two features $\in F$. Despite its simplicity, degree centrality can be very illuminating and can be considered a crude measure of whether a node is influential or not in the TMFG (i.e. variables mostly contributing to the system’s information flow). Once obtained $\deg(v)\;\forall v\in\textrm{TMFG}$, we rank these values in a descending order and we take the top $k$ central nodes, where $k$ is the cardinality of the features’ subset we want to consider.

To prove the effectiveness of the proposed methodology, we test the TFS algorithm against the current state-of-the-art counterpart in unsupervised, graph-based filter feature selection techniques, i.e. Infinite Feature Selection ($\textrm{Inf-FS}_{U}$) (Roffo et al., 2020). ¹¹1The Python implementation of $\textrm{Inf-FS}_{U}$ algorithm used in this paper can be reached at https://github.com/fullyz/Infinite-Feature-Selection. $\textrm{Inf-FS}_{U}$ represents features as nodes of a graph and relationships among them as weighted edges. Weights are computed as per in Equation 3. Each path of a given length over the network is seen as a potential set of relevant features. Therefore, varying paths and letting them tend to an infinite number permits the investigation of the importance of each possible subset of features. Based on this, assigning a score to each feature and ranking them in descendant order allows us to perform feature selection effectively. It is worth noting that $\textrm{Inf-FS}_{U}$ has a computational complexity equal to $\approx\mathcal{O}(n^{3})$. In contrast, TFS has a computational complexity equal to $\approx\mathcal{O}(n^{2})$, with $n$ being the number of the features.

Once defined the training, validation and test set for each benchmark dataset, model hyper-parameters (see Table 1) are optimised adopting a parallel grid search approach. For $\textrm{Inf-FS}_{U}$, the first hyper-parameter to be optimised is $\alpha$, which can take values between $0$ and $1$. To tune this parameter, we use a range of equally spaced realisations between $0.1$ and $1.0$, all at a distance of $0.1$. The second hyper-parameter to be optimised is $\theta$ and represents a regularisation factor which, in the original paper (Roffo et al., 2015), has a fixed value equal to $0.9$. Here, instead, we tune this parameter in the same way as $\alpha$. In the case of TFS, the first hyper-parameter to be optimised is the metric used in the building process of the initial fully connected graph. As reported in Section 2.3, we test three different metrics: (i) the Pearson’s estimator of the correlation coefficient, (ii) the Spearman’s rank correlation coefficient and (iii) the Energy coefficient. The second hyper-parameter to be tuned is a boolean value which regulates the chance to use the coefficients mentioned above in a squared form. It is worth mentioning that, if the Energy metric is chosen, the corresponding $\hat{\textbf{C}}$ is never squared since it already contains only positive values. The last hyper-parameter to be optimised is $\alpha$, and it should be considered only when the Energy coefficient is chosen as a metric. The meaning of this hyper-parameter is the same as its homologous in the $\textrm{Inf-FS}_{U}$ model. All the models are finally evaluated on feature subsets with cardinalities $\in[10,50,100,150,200]$.

For each hyper-parameter combination, a stratified $k$-fold cross-validation with $k=3$ is performed on the training set. The value of $k$ is chosen to take into account labels’ distributions reported in Appendix B. Results’ reproducibility and a fair comparison between models are guaranteed by fixing the random seed for each step of the training/validation/test pipeline.

The meaningfulness of each subset of features chosen by the two algorithms is evaluated based on the classification performance achieved by three classification algorithms: (i) Linear Support Vector Classifier (LinearSVM); (ii) k-Nearest Neighbors Classifier (KNN); (iii) Decision Tree Classifier (Decision Tree). LinearSVM is a sparse kernel-based method designed to convert non-linearly separable problems in the low-dimensional space, into linearly separable problems in the higher-dimensional space, thereby achieving classification (Han et al., 2022). KNN is a lazy learning algorithm, which classifies test instances evaluating their distance from the nearest k training samples stored in an n-dimensional space (where n is the number of dataset’s covariates) (Han et al., 2022). Finally, Decision Tree is a flowchart-like tree structure computed on training instances which classifies test samples tracing a path from the root to a leaf node holding the class prediction (Han et al., 2022). The inherently different nature of the three classifiers prevents from obtaining biased results for the two feature selection approaches. More details about chosen classifiers and their implementations are reported in Appendix D. Learning performances are evaluated based on three different metrics: (i) the Balanced Accuracy score (BA) (Mosley, 2013; Kelleher et al., 2020); (ii) the F1 score (F1); (iii) the Matthews Correlation Coefficient (MCC) (Baldi et al., 2000; Gorodkin, 2004; Jurman et al., 2012). We use the BA score for the hyper-parameters optimization process and as the reference metric to present results in Section 3. The BA score for the multi-class case is defined as:

TP is the number of outcomes where the model correctly classifies a sample as belonging to a positive class, when in fact it does belong to that class. TN is the number of outcomes where the model correctly classifies a sample as belonging to a negative class, when in fact it does not belong to that class. FP is the number of outcomes where the model incorrectly classifies a sample as belonging to a positive class, when in fact it does not belong to that class. FN is the number of outcomes where the model incorrectly classifies a sample as belonging to a negative class, when in fact it belongs to a positive class. $|Z|$ indicates the cardinality of the set of different classes.

General formulations for the F1 score, and the MCC are reported in Appendix F together with an extended version of results described later in this Section.

For each model and for each subset’s cardinality, the hyper-parameters configuration which maximises the BA score while minimising the number of parameters is applied to test datasets. In order to assess the robustness of our findings, we test if the results achieved by the two feature selection approaches are statistically different. To achieve this goal we use an improved version of the classic $5\times 2$ cv paired t-test (Dietterich, 1998). The test is constructed as follows. Given two classifiers $A$ and $B$ and a dataset $D$, $D$ is first randomly split into two balanced subsets $D_{1}$, $D_{2}$ (one for training and one for test). Both $A$ and $B$ are then estimated on $D_{1}$ and evaluated on $D_{2}$ obtaining performance measures $a_{1},b_{1}$. The roles of the datasets are then switched by estimating $A$ and $B$ on $D_{2}$ and evaluating on $D_{1}$ which results in further performance measures $a_{2},b_{2}$. The random division of $D$ is performed for a total of $5$ times, obtaining the matched performance evaluations $\{a_{1},b_{1}\},\{a_{2},b_{2}\}\dots,\{a_{10},b_{10}\}$. The test statistic t is then computed as follows:

where $d_{h}=a_{h}-b_{h}$ for $h=1,\dots,10$, is the difference between the matched performances metrics of the two classifiers, $\bar{d}=\frac{1}{10}\sum_{h=1}^{10}d_{h}$ and $\hat{\sigma}^{2}=\frac{1}{10}\sum_{h=1}^{10}(d_{h}-\bar{d})^{2}$. t follows a t-distribution with $9$ degrees of freedom and the null hypothesis is that the two classifiers $A$ and $B$ are not statistically different in their performances. Starting from this basic formulation of the $5\times 2$ cv paired t-test, we simply increased the number of iterations, making it a $15\times 2$ cv paired t-test, in order to increase the statistical robustness of the achieved results. Also in this case, results reproducibility is guaranteed through a strict control of random seeds.