IMMIGRATE: A Margin-Based Feature Selection Method with Interaction Terms

Traditional hypothesis-margin researches focus on obtaining large margins and feature selection. In this work, we show that the robustness of margins is also critical and can be measured using entropy. In addition, our approach provides clear mathematical formulations and explanations to uncover feature interactions, which is often lack in large hypothesis-margin based approaches. We design an algorithm, termed IMMIGRATE (Iterative max-min entropy margin-maximization with interaction terms), for training the weights associated with the interaction terms. IMMIGRATE simultaneously utilizes both local and global information and can be used as a base learner in Boosting. We evaluate IMMIGRATE in a wide range of tasks, in which it demonstrates exceptional robustness and achieves the state-of-the-art results with high interpretability.


Introduction
Feature selection is one of the most fundamental problems in machine learning and pattern recognition [1]. The Relief algorithm by Kira and Rendell [2] is one of the most successful feature selection algorithms. It can be interpreted as an online learning algorithm that solves a convex optimization problem with a hypothesis-margin-based cost function. Instead of deploying exhaustive or heuristic combinatorial searches, Relief decomposes a complex, global and nonlinear classification task into a simple and local one. Following the large hypothesis-margin principle for classification, Relief calculates the weights of features, which can be used for feature selection. Considering the binary classification in a set of samples P with two kinds of labels, the hypothesis-margin of an instance x is later formally defined in Gilad-Bachrach et al. [3] as 1 2 ( x − NM( x) − x − NH( x) ), where NH( x) denotes the "nearest hit," i.e., the nearest sample to x with the same label, while NM( x) denotes the "nearest miss", the nearest sample to x with the different label. The large hypothesis-margin principle has motivated several successful extensions of the Relief algorithm. For example, ReliefF [4] uses multiple nearest neighbors. Simba [3] recalculates the nearest neighbors every time the feature weights are updated. Yang et al. [5] consider global information to improve Simba. I-RELIEF [6] identifies the nearest hits and misses in a probabilistic manner, which forms a variation of hypothesis-margin. LFE [7] extends Relief from feature selection to feature extraction using local information. IM4E is proposed by Bei and Hong [8] to balance margin-quantity maximization and margin-quality maximization. Both approaches in Sun and Wu [7], Bei and Hong [8] use a variation of hypothesis-margin proposed in Sun and Li [6].
The Relief-based algorithms indirectly consider feature interactions by normalizing the feature weights [9], which, however, cannot directly reflect natural effects of associations and hence results in poor understanding on how feature interacts. For example, Relief and many of its extensions cannot tell whether a high weight of a certain feature is caused by its linear effect or its interaction with other features [9]. Furthermore, these methods cannot directly reveal and measure the impact of the interaction terms on classification results.
To this end, we propose the Iterative Max-MIn entropy marGin-maximization with inteRAction TErms algorithm (IMMIGRATE, henceforth). IMMIGRATE directly measures the influence of feature interactions and has the following characteristics. First, when defining hypothesis-margin, we introduce a new trainable quadratic-Manhattan measurement to capture interaction terms, which measures the interaction importance directly. Second, we take advantage of the margin stability by measuring the underlying entropy based on the distribution of instances. Third, we derive an iterative optimization algorithm to efficiently minimize the cost function. Fourth, we design a novel classification method that utilizes the learned quadratic-Manhattan measurement to predict the class of a new instance. Fifth, we design a more powerful approach (i.e., Boosted IMMIGRATE) by using IMMIGRATE as the base learner of Boosting [10]. Sixth, to make IMMIGRATE efficient for analyzing high-dimensional datasets, we take advantage of IM4E [8] to obtain an effective initialization.
The rest of the paper is organized as follows. Section 2 explains the foundation of the Relief algorithm, and Section 3 introduces the IMMIGRATE algorithm. Section 4 summarizes and discusses our experiments with different datasets, showing that IMMIGRATE achieves the state-of-the-art results, and Boosted IMMIGRATE outperforms other boosting classifiers significantly. The computation time of IMMIGRATE is comparable to other popular feature selection methods that consider interaction terms. Section 5 concludes the article with comparisons with related works and a short discussion.

Review: the Relief Algorithm
We first introduce a few notations used throughout the paper: x i ∈ R A as the i-th instance in the training set P; y i as the class label of x i ; N as the size of P; A as the number of features (i.e., attributes); w as the feature weight vector; and | x i | as a vector where absolute value operation is element-wise. Relief [2] iteratively calculates the feature weights in w (Algorithm 1). The higher a feature weight is, the more relevant the corresponding feature is. After the calculation of feature weights, a threshold is chosen to select relevant features. Relief can be viewed as a convex optimization problem that minimizes the cost function in Equation 1: where M( N) is a user defined number of randomly chosen training samples, NH( x) is the nearest "hit" (from the same class) of x; NM( x) is the nearest "miss" (from a different class) of x; and w T x n − NH( x n ) is the weighted Manhattan distance. Denote u = ∑ M n=1 x n − NH( x n ) − x n − NM( x n ) . Minimizing the cost function of Relief (1) can be solved using the Lagrange multiplier method and the Karush-Kuhn-Tucker conditions [11] to get a closed-form solution: w = (− u) + / (− u) + 2 , where ( a) + truncates the negative elements to 0. This solution to the original Relief algorithm is important for understanding the Relief-based algorithms.

Algorithm 1 The Original Relief Algorithm
N: the number of training instances.
A: the number of features(i.e. attributes).
M: the number of randomly chosen training samples to update feature weight w.
Randomly select an instance x i and find its NH( x i ) and NM( x i ). Update the feature weights by where the square operation is element-wise. Return: w.

IMMIGRATE Algorithm
Without loss of generality, we establish the IMMIGRATE algorithm in a general binary classification setting. This formulation can be easily extended to handle multi-class classification problems. Let the whole data set be P = {z n | z n = ( x n , y n ), x n ∈ R A , y n = ±1} N n=1 ; the hit index set of x n be H n = {j | z j ∈ P, y j = y n & j = n}, and the miss index set of x n be M n = {j | z j ∈ P, y j = y n }.

Hypothesis-Margin
Given a distance d( x i , x j ) between two instances, x i and x j , a hypothesis-margin [3] is defined as , where x h and x m represent the nearest hit and nearest miss for instance x n , respectively. We adopt the probabilistic hypothesis-margin defined by Sun and Li [6] as where α n,h ≥ 0, β n,m ≥ 0, ∑ h∈H n α n,h = 1, ∑ m∈M n β n,m = 1, for ∀ n ∈ {1, · · · , N}. As in the above design, the hidden random variable α n,h represents the probability that x h is the nearest hit of instance x n , while β n,m indicates the probability that x m is the nearest miss of instance x n . In the rest of the paper, for conciseness, we will use margin to indicate hypothesis-margin.

Entropy to Measure Margin Stability
The distributions of hits and misses can be used to evaluate the stability of margins (i.e., margin quality). A more stable margin can be obtained by considering the distributions of instances with the same or different labels with respect to the target instance. A margin is deemed stable if it will not be greatly reduced by changes to only a few neighbors of the target instance. Considering an instance x n , its probabilities {α n,h } and {β n,m } represent the distributions of its hits and misses, respectively. We can use the hit entropy E hit ( x n ) = − ∑ h∈H n α n,h log α n,h and miss entropy E miss ( x n ) = − ∑ m∈M n β n,m log β n,m to evaluate the stability of x n 's margin. The following two scenarios help explain the intuition of using these entropy. Scenario A: all neighbors are distributed evenly around the target instance; scenario B: the neighbor distribution is highly uneven. An extreme example for scenario B is that one instance is quite close to the target and the rest are quite far away from the target. An easy experiment to test the stability is to discard one instance from the system and to check how it influences the margin. In scenario A, if the closest neighbor (no matter if it is hit or miss) is discarded, the margin changes only slightly because there are many other hits/misses evenly distributed around the target. In scenario B, if the closest neighbor is a miss, its removal can increase the margin significantly. On the contrary, if the closest neighbor is a hit, removing it can decrease the margin significantly. Intuitively speaking, hits prefer scenario A and misses favor scenario B.
Since scenarios A and B correspond to high and low entropy, respectively, the margin can benefit from a large hit entropy E hit (e.g., scenario A) and a low miss entropy E miss (e.g., scenario B). We can set up a framework to maximize the hit entropy and minimize the miss entropy, which is equivalent to make the margin in Equation 2 the most stable. Bei and Hong [8] use the term max-min entropy principle to describe the process that maximizes the hit entropy and minimize the loss entropy to maximize the margin quality. The process of stabilizing margin is an extension of the large margin principle.

Quadratic-Manhattan Measurement
We extend the margin in Equation 2 by using a new quadratic-Manhattan measurement defined as: where W is a non-negative symmetric matrix (element-wise non-negative) with its Frobenius norm W F = 1. The quadratic-Manhattan measurement is a natural extension of the weight vector, and the distance defined in Equation 3 is a natural extension of the weighted Manhattan distance in Equation 1. Off-diagonal elements in W capture feature interactions and diagonal elements in W capture main effects. To understand why quadratic-Manhattan measurement can capture the influence of interactions, we observe that the effect of element w a,b (a = b) in W enters into (3) as the coefficient for the combination of the a-th and b-th elements in vector x i − x j . In Relief-based algorithms, the weighted Manhattan distance Equation 1 can be equivalently captured by the feature weight update equation in Algorithm 1. Similarly, w a,b can be updated using the combination of the a-th and b-th features based on a randomly given instance. We thus define our new margin using the quadratic-Manhattan measurement as

IMMIGRATE
We design the following cost function to maximize our new margin, and simultaneously, the hit entropy and miss entropy are optimized.
where E miss (z n ) = − ∑ m∈M n β n,m log β n,m , E hit (z n ) = − ∑ h∈H n α n,h log α n,h , and σ is a hyperparameter that can be tuned via internal cross-validation. We also design the following optimization procedure containing two iterative steps to find W that minimizes the cost function. The framework starts from a randomly initialized W and stops when the change of cost function is less than a preset limit or the iteration number reaches a preset threshold. In practice, we find that it typically takes less than 10 iterations to stop and obtain good results. Based on our experiments, different initialization of W does not influence the results of the iterative optimization. The computation time of IMMIGRATE is comparable to other interaction related methods such as SODA [12], hierNet [13].
As depicted by the flow-chart in Figure 1, the IMMIGRATE algorithm iteratively optimizes the cost function Equation 5. It starts with a random initiation satisfying certain boundary conditions, and proceeds to iterate the two steps as detailed below in Algorithm 2.
is small enough or the iteration indicator t reaches a preset limit.
The Hessian matrix of C w.r.t. probability pair (α n,h , β n,m ) is: Since α n,h , β n,m > 0, the determinant of the Hessian matrix is negative, where a saddle point is found in the (α n,h , β n,m ) space. Therefore, the cost function C achieves its local minimum and local maximum w.r.t. α n,h and β n,m , respectively.
Let the ψ i 's and µ i 's be the eigenvectors and eigenvalues of Σ, respectively, so that where Proof. Since W is a distance metric matrix, it is symmetric and positive-semidefinite. Let λ 1 ≥ λ 2 ≥ · · · ≥ λ A ≥ 0 be eigenvalues of W, then the eigen-decomposition of W is where P is an orthogonal matrix, and The constraint W 2 F = 1 can be simplified as: Let us rearrange Equation 5 as: Then, Equation 5 can be further simplified as: where Σ = ∑ N n=1 Σ n,H − Σ n,M and Σ n,H = ∑ h∈H n α n,h x n − x h x n − x h T , Σ n,M = ∑ m∈M n β n,m x n − x m x n − x m T . The orthogonality condition can be ignored because this condition is required in the constraint. The Lagrangian for the optimization problem in Equation 12 is easy to obtain: Differentiating L with respect to φ i yields: Denote where µ i = −2λ φ i 2 2 . Thus, ψ i and µ i are an eigenvector and eigenvalue of Σ, respectively. Equation 12 can be simplified to be Note that Equation 16 is exactly the same as the original Relief Algorithm (Algorithm 1): where ( a) + = [max(a 1 , 0), max(a 2 , 0), · · · , max(a I , 0)], and φ i = √ η i ψ i . It is also easy to see that the updated W is a distance metric.

Weight Pruning
Some previous Relief-based algorithms offer options to remove weights lower than a preset threshold. IMMIGRATE offers a similar option to prune small weights: set small elements in W to 0. By default, we use a threshold to prune small weights to 0, where W should be normalized w.r.t. Frobenius norm after the pruning.

Predict New Samples
A prediction rule based on the learned weight matrix W can be formulated as: where z = ( x , y ) is a new instance, c denotes the class andŷ is the predicted label. This prediction method assigns a new instance to a class that maximizes its hypothesis-margin using the learned weight matrix W, which makes it more stable than the k-NN method used in the traditional Relief-based algorithms.

IMMIGRATE in Ensemble Learning
Boosting [10,14,15] has been widely used to create ensemble learners that produce the state-of-the-art results in many tasks. Boosting combines a set of relatively weak base learners to create a much stronger learner. To use IMMIGRATE as the base classifier in the AdaBoost algorithm [14], we modify the cost function Equation 5 to include sample weights and use the modified version in the boosting iterations. We name the algorithm BIM, standing for Boosted IMMIGRATE (Refer to Equation 19 and Algorithm 3 for more details about BIM. BIM schedules the adjustment of the hyperparameter σ in its boosting iterations. It starts with σ being a predefined σ max and gradually reduces σ by multiplying it with (σ min /σ max ) 1/T at each interaction until reaching σ min , where T is a predefined maximum number of boosting iterations.
where E miss (z n ) = − ∑ m∈M n β n,m log β n,m , E hit (z n ) = − ∑ h∈H n α n,h log α n,h , ∑ N n=1 D( x n ) = 1, and D( x n ) ≥ 0, ∀ n Algorithm 3 The BIM Algorithm T: the number of classifiers for BIM. Input : a training dataset {z n = ( x n , y n )} n=1,··· ,N . Initialization : for each x n , set D 1 ( x n ) = 1/N. for t := 1 to T do Limit max number of iteration of IMMIGRATE less than preset. Train weak IMMIGRATE classifier h t (x) using a chosen σ t and weights D t (x) by Equation 19.
Compute the error rate t as Discard h t , T = T − 1 and continue .

IMMIGRATE for High-Dimensional Data Space
When applied to high-dimensional data, IMMIGRATE can incur a high computational cost because it considers the interactions between every feature pair. To reduce the computational cost, we first use IM4E [8] to learn a feature weight vector, which is used to initialize the diagonal elements of W in the proposed quadratic-Manhattan measurement. We also use the learned feature weight vector to help pre-screen the features, and keep only those with weights above a preset limit. In the remaining computation, we only model interactions between those chosen features. The discarded features after pre-screening can be added back empirically based on the need of a specific application. We term this procedure IM4E-IMMIGRATE, which is effective and computationally efficient. It can also be boosted (Boosted IM4E-IMMIGRATE) to be stronger.

Experiments
In our experiments, all continuous features are normalized with mean zero and unit variance.
And cross-validation is used here to compare the performances of various approaches. We have implemented IMMIGRATE in R and MATLAB. The R package is available at https://CRAN.R-project.org/package=Immigrate, and the MATLAB version is available at https://github.com/RuzhangZhao/Immigrate-MATLAB-. Both IMMIGRATE and BIM can be accelerated by parallel computing as their computations are matrix-based.

Synthetic Dataset
We first test the robustness of the IMMIGRATE algorithm using a synthesized dataset where we have two interacting features following Gaussian distributions in a binary classification setting. The simulated dataset contains 100 samples from one class governed by a Gaussian distribution with mean (4, 2) T and the covariance matrix 1 0.5 0.5 1 and another 100 samples from the other class governed by a Gaussian distribution with mean (6, 0) T and the same covariance matrix. In addition, we add noises following a Gaussian distribution with mean (8, −2) T and the covariance matrix 8 4 4 8 to the fist class, and add noises following a Gaussian distribution with mean (2, 4) T and the same covariance matrix to the second class. Figure 2 shows a scatter plot of the synthesized dataset containing 10% samples from the noise distributions. The slope of the orange dotted line in Figure 2 is 1, which separates data with different labels. The noises are included to disturb the detection of the interaction term. The noise level starts from 5%, and gradually increases by 5% to 50%. As the baseline, we apply logistic regression and observe that the t-test p-value of the interaction coefficient increases from 3 × 10 −11 to 7 × 10 −5 and 0.7 when the noise level increases from 0% to 10% and 50%. Local Feature Extraction (LFE, Sun and Wu [7]) is a Relief-based algorithm which considers interaction terms indirectly, though the interaction information is only used for feature extraction. We run IMMIGRATE and LFE on the synthesized datasets and compare the weights of the interaction term between features 1 and 2 in Figure 3, which shows IMMIGRATE is more robust than LFE.
Whenever possible, we use the settings of the aforementioned methods reported in their original papers: LMNN uses 3-NN classifier; Relief and Simba use Euclidean distance and 1-NN classifier; ReliefF uses Manhattan distance and k-NN classifier (k=1,3,5 is decided by internal cross-validation); in SODA, gam (=0, 0.5, 1) is determined by internal cross-validation and logistic regression is used for prediction. The IM4E algorithm has two hyperparameters λ and σ. We fix λ = 1 as it has no actual contribution and tune σ as suggested by Bei and Hong [8]. Hence, the IMMIGRATE algorithm only has one hyperparameter σ. When tuning σ, we gradually decrease σ from σ 0 = 4 by half each time until it is not larger than 0.2. The preset limit for weight pruning is 1/A, where A is the number of features. Furthermore, the preset iteration number is chosen to be 10. For each dataset, σ and whether weight pruning is applied are determined by the best internal cross-validation results. For BIM, we use σ max = 4, σ min = 0.2, and the maximal number of boosting iterations T is 100. The preset threshold in IM4E-IMMIGRATE is 2/A.
We repeat ten-fold cross-validation ten times for each algorithm on each dataset, i.e., 100 trials are carried out. When comparing two algorithms (i.e., A vs. B), we calculate the paired Student's t-test using the results of 100 trials. First, the null hypothesis is there is no difference between the performances of A and those of B. When the p-value is larger than the significant level cutoff 0.05, we say A "Tie" B, which means there is no significant difference between their performances. When the p-value is smaller than the significant level cutoff 0.05, the second null hypothesis is the performances of B are no worse than those of A. When the new p-value is smaller than the significant level cutoff 0.05, we say A "wins", which means A on average performs significantly better than B on this dataset, and vice versa.

Gene Expression Datasets
Gene expression datasets typically have thousands of features. We use the following five gene expression datasets for feature selections: GLI [26], Colon (COL) [27], Myeloma (ELO) [28], Breast (BRE) [29], Prostate (PRO) [30]. All datasets have more than 10,000 features. Refer to Table A1 in Appendix A for details of all datasets.

UCI Datasets
We also carry out an extensive comparison using many UCI datasets [31]:  Table A1 for the full names and links for those datasets. If a dataset has more than two classes, we use two classes with the largest sample size. In addition, we use three large-scale datasets: CRO * , ELE * , WAV * .
We perform ten-fold cross-validation ten times. Tables 2 for IMMIGRATE and Table 3 for BIM show the average accuracies on the corresponding datasets. In Table 2, the last row "(W,T,L)" indicates the number of times IMMIGRATE (IGT) and BIM W,T,L (win,tie,loss) when compared with each algorithm separately by using the paired Student's t-test with the significance level of α = 0.05. The comparison results are also summarized in Figure 4 (bottom subplot), where the first 17 items (black) indicate the results for IMMIGRATE while the last three items (blue) indicate the results for BIM.
Although IMMIGRATE or BIM is not always the best, they outperform other methods significantly in one-to-one comparisons in terms of cross-validation results. Figure 4 (bottom subplot, black part) and Table 2 show that IMMIGRATE achieves the state-of-the-art performance as the base classifier while Figure 4 (bottom subplot, blue part) and Table 3 show BIM achieves the state-of-the-art performance as the boosted version. To visualize the feature selection results of our approaches, we plot the feature weight heat maps of four datasets (GLA, LYM, SMR and STA) in Appendix B Figure A1.
the number of times that the Boosted IMMIGRATE (BIM) wins/ties/losses an existing algorithm according to the paired t-test on the cross-validation results.

Related Works
In many recent publications, Relief-based algorithms and feature selection with interaction terms have been well explored. Some methods are reviewed here to show the connection and differences with our approach. The hypothesis-margin definition in Equation 2 adopted in this work is also used in some previous studies, such as Bei and Hong [8]. However, Bei and Hong [8] do not consider the interactions between features. Our work provides a measurable way to show the influence of each feature interaction.
Sun and Wu [7] propose local feature extraction (LFE) method, which learns linear combination of features for feature extraction. LFE explores the information of feature interaction terms indirectly, which is partly our aim. However, LFE does not consider global information or margin stability, which results in significant differences in the cost function and the optimization procedures.
Our quadratic-Manhattan measurement defined in Equation 3 is related to the Mahalanobis metric used in previous works on metric learning, such as Large Margin Nearest Neighbor (LMNN) [21]. Weinberger and Saul [21] use semi-definite programming for learning distance metric in LMNN. LMNN and our approach are both based on K-Nearest Neighbors. A major difference is that our quadratic-Manhattan measurement has matrix W to be non-negative and symmetric (element-wise non-negative) with its Frobenius norm W F = 1, whereas metric learning only requires its matrix to be symmetric semi-positive definite. Actually, the non-negative element requirement of W provides IMMIGRATE a high intepretability, where items in matrix indicate interaction importance.
Quadratic-Manhattan measurement serves well in the classification task and offers a direct explanation about how features, in particular, feature interaction terms, contribute to the classification results.

Conclusions and Discussion
In this paper, we propose a new quadratic-Manhattan measurement to extend the hypothesis-margin framework, based on which a feature selection algorithm IMMIGRATE is developed for detecting and weighting interaction terms. We also develop its extended versions, Boosted IMMIGRATE (BIM) and IM4E-IMMIGRATE. IMMIGRATE and its variants follow the principle of maximizing stable hypothesis-margin and are implemented via a computationally efficient iterative optimization procedure. Extensive experiments show that IMMIGRATE outperforms state-of-the-art methods significantly, and its boosted version BIM outperforms other boosting-based approaches. In conclusion, compared with other Relief-based algorithms, IMMIGRATE mainly has the following advantages: (1) both local and global information are considered; (2) interaction terms are used; (3) robust and less prone to noise; (4) easily boosted. The computation time of IMMIGRATE variants is comparable to other methods able to detect interaction terms.
There are some limitations for IMMIGRATE and we discuss some directions of improving the algorithm accordingly. First, in Section 3.4.3, small weights are removed to obtain sparse solutions using some cutoffs directly, which is hard to do inference for the obtained weights. Penalty terms such as the l 1 -or l 2 -penalty are usually applied to shrink and select important weights. We suggest that our cost function Equation 5 can be modified to include such a penalty term to replace the process of weight pruning in Section 3.4.3. Second, although IMMIGRATE is efficient, it still costs much time to compute data with large size. To further improve the computational efficiency of IMMIGRATE for large-scale datasets, we can improve training by using well selected prototypes [32], which, as a subset of the original data, are representative but with noisy and redundant samples removed. Third, IMMIGRATE only considers pair-wise interactions between features. Interactions among multiple features can play important roles in real applications, [33,34]. Our work provides a basis for developing new algorithms to detect multi-feature interactions. For example, people can use tensor form to consider weights for multi-feature interactions. Fourth, although our iterative optimization procedure is efficient, it achieves ad hoc solutions with no guarantee of reaching the global optimum. It remains an open challenge to develop better optimization algorithms. Finally, the selection of an appropriate σ currently relies on internal cross-validation, which cannot uncover the underlying properties of σ. A better strategy may be developed by rigorously investigating the theoretical contributions of σ.