Efficient Approximate Solutions to Mutual Information Based Global Feature Selection

TL;DR

This paper introduces two novel global approximation algorithms, TPower and LowRank, for efficiently solving the NP-hard problem of mutual information-based feature selection under conditional independence assumptions, demonstrating superior performance.

Contribution

The paper presents two new algorithms, TPower and LowRank, that provide effective approximations for the intractable CMI-based feature selection problem, improving over existing methods.

Findings

TPower and LowRank algorithms achieve high-quality feature selection.

Experimental results show superior performance over existing methods.

Algorithms are effective across multiple datasets.

Abstract

Mutual Information (MI) is often used for feature selection when developing classifier models. Estimating the MI for a subset of features is often intractable. We demonstrate, that under the assumptions of conditional independence, MI between a subset of features can be expressed as the Conditional Mutual Information (CMI) between pairs of features. But selecting features with the highest CMI turns out to be a hard combinatorial problem. In this work, we have applied two unique global methods, Truncated Power Method (TPower) and Low Rank Bilinear Approximation (LowRank), to solve the feature selection problem. These algorithms provide very good approximations to the NP-hard CMI based feature selection problem. We experimentally demonstrate the effectiveness of these procedures across multiple datasets and compare them with existing MI based global and iterative feature selection procedures.

Tables

Table 1. TABLE I: Time complexities for the Global approximate solutions for BQP in number of features n 𝑛 n

Linear	Spectral	SDP	LowRank*	TPower**
$O\left(nk\right)$	$O\left(n^2\right)$	$O\left(n^{4.5}\right)$	$O\left(n^{\left(d+1\right)}\right)$	$O\left(t{n}^2\right)$
*$d$ is approximation rank,  **$t$ is number of iterations

Table 2. TABLE II: Datasets details: n 𝑛 n is number of features, m 𝑚 m is number of samples, c 𝑐 c is number of categories, Error : is average cross validation error (%) using all features.

Data	$n$	$m$	$c$	Error	Ref.
Arrhythmia	258	420	2	31.1	[20]
Colon	2000	62	2	37.0	[21]
Gisette	4995	6000	2	2.5	[20]
Leukemia	7070	72	2	26.4	[21]
Lung	325	73	7	9.6	[21]
Lymphoma	4026	96	9	81.3	[21]
Madelon	500	2000	2	45.5	[20]
Multi-Feat	649	2000	10	1.5	[20]
Musk2	166	6598	2	4.6	[20]
OptDigits	62	3823	10	3.3	[20]
Promoter	57	106	2	26.0	[20]
Spambase	57	4601	2	7.5	[20]
Waveform	21	5000	3	13.1	[20]

Table 3. TABLE III: Comparison of TPower with other algorithms. The table values measure the difference in average classification accuracies of TPower with other algorithms. w 𝑤 w , t 𝑡 t and l 𝑙 l indicate one-sided paired t-test results. The last row displays the total number of Wins( W 𝑊 W ), Ties( T 𝑇 T ) and Loss( L 𝐿 L ). N/A indicates comparison data was unavailable for large datasets using SDP .

Data	MaxRel	MRMR	JMI	QPFS	Spectral	SDP	LowRank
Arrythmia	-0.37 $\pm$ 1.4 t	0.32 $\pm$ 1.0 l	0.02 $\pm$ 1.0 l	0.20 $\pm$ 1.8 l	-0.08 $\pm$ 1.1 t	-0.18 $\pm$ 1.0 w	-0.05 $\pm$ 0.8 l
Colon	-7.28 $\pm$ 4.6 w	-4.42 $\pm$ 4.2 w	-2.47 $\pm$ 3.8 w	-6.70 $\pm$ 4.6 w	-0.60 $\pm$ 2.8 w	N/A	4.03 $\pm$ 4.5 l
Gisette	-1.32 $\pm$ 0.6 w	0.00 $\pm$ 0.7 w	-1.12 $\pm$ 0.6 w	-1.38 $\pm$ 0.7 w	-1.26 $\pm$ 0.6 w	N/A	0.33 $\pm$ 0.6 l
Leukemia	0.11 $\pm$ 1.4 w	1.40 $\pm$ 1.6 l	1.59 $\pm$ 1.8 l	0.41 $\pm$ 1.1 t	-0.03 $\pm$ 0.6 w	N/A	1.49 $\pm$ 1.3 l
Lung	-9.43 $\pm$ 4.1 w	-2.52 $\pm$ 4.2 w	-3.83 $\pm$ 4.2 w	0.60 $\pm$ 2.8 l	-0.88 $\pm$ 2.2 w	-0.88 $\pm$ 2.1 w	-1.59 $\pm$ 2.4 w
Lymphoma	-2.76 $\pm$ 4.8 w	3.35 $\pm$ 4.7 l	2.93 $\pm$ 5.3 l	4.99 $\pm$ 3.3 l	-1.86 $\pm$ 2.5 w	N/A	3.29 $\pm$ 4.2 l
Madelon	0.32 $\pm$ 0.5 l	0.80 $\pm$ 0.9 l	0.01 $\pm$ 0.4 w	-0.22 $\pm$ 0.7 w	-0.01 $\pm$ 0.6 w	0.15 $\pm$ 0.6 l	-0.11 $\pm$ 0.4 t
MultiFeatures	0.02 $\pm$ 0.3 w	0.24 $\pm$ 0.3 l	0.17 $\pm$ 0.3 l	-0.42 $\pm$ 0.3 w	0.10 $\pm$ 0.3 w	0.11 $\pm$ 0.3 w	0.01 $\pm$ 0.3 l
Musk2	-0.45 $\pm$ 0.6 w	-0.22 $\pm$ 0.7 w	-0.18 $\pm$ 0.5 w	-0.31 $\pm$ 0.6 w	0.06 $\pm$ 0.4 w	0.03 $\pm$ 0.5 w	0.05 $\pm$ 0.4 w
OptDigits	-0.19 $\pm$ 0.5 w	-0.01 $\pm$ 0.6 t	0.16 $\pm$ 0.6 l	-0.65 $\pm$ 1.0 w	0.03 $\pm$ 0.3 l	-2.53 $\pm$ 13.0 w	0.08 $\pm$ 0.4 l
Promoter	0.73 $\pm$ 3.0 l	-0.04 $\pm$ 3.2 w	0.19 $\pm$ 3.0 l	-1.29 $\pm$ 3.8 w	-0.48 $\pm$ 2.8 w	-0.56 $\pm$ 2.9 w	-0.17 $\pm$ 3.1 w
Spambase	-0.34 $\pm$ 0.3 w	0.06 $\pm$ 0.2 l	-0.23 $\pm$ 0.3 w	0.03 $\pm$ 0.4 l	-0.09 $\pm$ 0.3 w	-0.10 $\pm$ 0.3 w	0.02 $\pm$ 0.1 l
Waveform	-0.13 $\pm$ 0.3 w	0.06 $\pm$ 0.3 l	-0.01 $\pm$ 0.0 t	0.04 $\pm$ 0.2 t	0.04 $\pm$ 0.1 t	0.00 $\pm$ 0.2 t	-0.01 $\pm$ 0.2 t
#$W$/$T$/$L$:	10/1/2	5/1/7	6/1/6	7/2/4	10/2/1	7/1/1	3/2/8

Table 4. TABLE IV: Comparison of LowRank with other algorithms. Table structure similar to Table III

Data	MaxRel	MRMR	JMI	QPFS	Spectral	SDP	TPower
Arrythmia	-0.32 $\pm$ 1.3 l	0.36 $\pm$ 1.0 l	0.07 $\pm$ 1.0 l	0.25 $\pm$ 1.7 l	-0.03 $\pm$ 1.1 t	-0.13 $\pm$ 1.0 w	0.05 $\pm$ 0.8 w
Colon	-11.31 $\pm$ 4.7 w	-8.45 $\pm$ 4.3 w	-6.50 $\pm$ 3.9 w	-10.73 $\pm$ 5.3 w	-4.63 $\pm$ 5.0 w	N/A	-4.03 $\pm$ 4.5 w
Gisette	-1.65 $\pm$ 0.5 w	-0.32 $\pm$ 0.5 w	-1.44 $\pm$ 0.7 w	-1.70 $\pm$ 0.6 w	-1.58 $\pm$ 0.6 w	N/A	-0.33 $\pm$ 0.6 w
Leukemia	-1.39 $\pm$ 1.4 w	-0.09 $\pm$ 1.5 t	0.09 $\pm$ 1.7 t	-1.09 $\pm$ 1.3 w	-1.52 $\pm$ 1.2 w	N/A	-1.49 $\pm$ 1.3 w
Lung	-7.83 $\pm$ 4.1 w	-0.92 $\pm$ 4.5 w	-2.23 $\pm$ 4.5 w	2.19 $\pm$ 3.7 l	0.71 $\pm$ 2.5 l	0.71 $\pm$ 2.3 l	1.59 $\pm$ 2.4 l
Lymphoma	-6.06 $\pm$ 3.5 w	0.06 $\pm$ 2.1 l	-0.36 $\pm$ 2.1 l	1.70 $\pm$ 2.5 l	-5.15 $\pm$ 3.8 w	N/A	-3.29 $\pm$ 4.2 w
Madelon	0.43 $\pm$ 0.5 l	0.91 $\pm$ 0.8 l	0.12 $\pm$ 0.5 l	-0.11 $\pm$ 0.7 w	0.10 $\pm$ 0.6 l	0.26 $\pm$ 0.5 l	0.11 $\pm$ 0.4 t
MultiFeatures	0.01 $\pm$ 0.4 l	0.23 $\pm$ 0.3 l	0.16 $\pm$ 0.4 l	-0.43 $\pm$ 0.3 w	0.10 $\pm$ 0.4 l	0.10 $\pm$ 0.4 l	-0.01 $\pm$ 0.3 w
Musk2	-0.50 $\pm$ 0.4 w	-0.27 $\pm$ 0.7 w	-0.23 $\pm$ 0.5 w	-0.36 $\pm$ 0.5 w	0.02 $\pm$ 0.3 l	-0.02 $\pm$ 0.4 l	-0.05 $\pm$ 0.4 l
OptDigits	-0.26 $\pm$ 0.6 w	-0.09 $\pm$ 0.4 w	0.08 $\pm$ 0.4 l	-0.72 $\pm$ 1.2 w	-0.04 $\pm$ 0.4 t	-2.61 $\pm$ 13.1 w	-0.08 $\pm$ 0.4 w
Promoter	0.90 $\pm$ 3.1 l	0.13 $\pm$ 2.7 w	0.35 $\pm$ 2.9 w	-1.13 $\pm$ 3.7 w	-0.31 $\pm$ 2.8 t	-0.40 $\pm$ 2.5 t	0.17 $\pm$ 3.1 l
Spambase	-0.36 $\pm$ 0.3 w	0.04 $\pm$ 0.2 l	-0.24 $\pm$ 0.3 w	0.02 $\pm$ 0.4 t	-0.11 $\pm$ 0.3 w	-0.12 $\pm$ 0.3 w	-0.02 $\pm$ 0.1 w
Waveform	-0.12 $\pm$ 0.4 w	0.07 $\pm$ 0.1 l	0.00 $\pm$ 0.2 t	0.05 $\pm$ 0.1 t	0.05 $\pm$ 0.1 t	0.02 $\pm$ 0.1 t	0.01 $\pm$ 0.2 t
#$W$/$T$/$L$:	9/0/4	6/1/6	6/2/5	8/2/3	5/4/4	3/2/4	8/2/3

Equations

$$p(X,Y)$$

$$=p(X_{\tilde{\mathbb{S}}})p(X_{\mathbb{S}}|X_{\tilde{\mathbb{S}}})p(Y|X_{\mathbb{S}})$$

$$I(X_{\mathbb{S}};Y)$$

$$=I(X_{{\mathbb{S}}_i};Y)+I(X_i;Y|X_{{\mathbb{S}}_i})$$

$$=I(X_{{\mathbb{S}}_i};Y)+I(X_i;Y)$$

$$-I(X_i;X_{{\mathbb{S}}_i})+I(X_i;X_{{\mathbb{S}}_i}|Y)$$

$$=I(X_{{\mathbb{S}}_i};Y)+I(X_i;Y)-H(X_{{\mathbb{S}}_i})+H(X_{{\mathbb{S}}_i}|X_i)$$

$$+H(X_{{\mathbb{S}}_i}|Y)-H(X_{{\mathbb{S}}_i}|X_i,Y)$$

$$=I(X_i;Y)+H(X_{{\mathbb{S}}_i}|X_i)-H(X_{{\mathbb{S}}_i}|X_i,Y)$$

$$I(X_i{X}_{{\mathbb{S}}_i};Y)$$

$$-\sum _{j\in {\mathbb{S}}_i}H(X_j|X_i,Y)$$

$$=I(X_i;Y)+\sum _{j\in {\mathbb{S}}_i}I(X_j;Y|X_i)$$

$$\mathbb{S}$$

$$\underset{\mathbf{x}}{\max }\{{\mathbf{x}}^{\top }\mathbf{Qx}\}\text{s.t.}\mathbf{x}\in \{0,1{\}}^n,||\mathbf{x}|{|}_1=k,$$

$$\underset{\mathbf{x}}{\max }$$

$$x_i,w_{ij}\in \{0,1\},\forall 1\le i,j\le n.$$

$$\text{Linear}=\underset{\mathbf{x}}{\max }$$

$$s.t.$$

$${\lambda }_{\max }(A,k)=\max {\mathbf{x}}^{\top }A\mathbf{x},s.t.\Vert \mathbf{x}\Vert =1,\Vert \mathbf{x}{\Vert }_0\le k$$

$${\text{BQP}}_b=$$

$$||\mathbf{x}|{|}_1=||\mathbf{y}|{|}_1=k$$

$$q:\{0,1{\}}^n\to (-\infty ,0],\text{where}$$

$$q(\mathbf{x})=\text{BQP}-||\mathbf{Q}|{|}_1$$

$$\underset{\mathbf{x}}{\mathrm{argmax}}\text{BQP}\equiv \underset{\mathbf{x}}{\mathrm{argmax}}q(\mathbf{x})$$

$$||\mathbf{Q}|{|}_1=Q^0+Q^1+Q^2\text{where},$$

$$Q^0\leftarrow \sum _{i,j|{\bar{x}}_i+{\bar{x}}_j=0}{\mathbf{Q}}_{ij}$$

$$Q^1\leftarrow \sum _{i,j|{\bar{x}}_i+{\bar{x}}_j=1}{\mathbf{Q}}_{ij}$$

$$Q^2\leftarrow \sum _{i,j|{\bar{x}}_i+{\bar{x}}_j=2}{\mathbf{Q}}_{ij}\equiv {\bar{\mathbf{x}}}^{\top }\mathbf{Q}\bar{\mathbf{x}}$$

$$||\mathbf{Q}|{|}_1$$

$$Q_{LP2}^{\ast }$$

$$=||\mathbf{Q}\bar{\mathbf{x}}|{|}_1$$

$$\ge ||\mathbf{Q}{\mathbf{x}}^{\ast }|{|}_1$$

$$\ge {\mathbf{x}}^{\ast \top }\mathbf{Q}{\mathbf{x}}^{\ast }$$

$$=Q^{\ast }$$

$$2q({\mathbf{x}}^{\ast })\le q(\bar{\mathbf{x}})$$

$${\mathbf{x}}^{\ast \top }\mathbf{Q}{\mathbf{x}}^{\ast }$$

$$2{\mathbf{x}}^{\ast \top }\mathbf{Q}{\mathbf{x}}^{\ast }$$

$$2q({\mathbf{x}}^{\ast })$$

Full text

I Introduction

High dimensional data can pose a significant challenge to learning methods due to the curse of dimensionality [1]. Feature selection is a prominent dimensionality reduction technique that selects a small subset of features based on certain relevancy criteria. Apart from reducing data dimensionality, feature selection provides insights into the data, prevents over-fitting and reduces computational costs for learning, which ultimately results in better learned models.

Depending on whether there is label information available, feature selection can be classified into two categories: supervised and unsupervised. Supervised feature selection procedures are broadly classified into three groups, wrapper, filter and embedded methods [2]. Wrapper procedures select features for a specific learning model. Filter methods on the other hand are classifier agnostic. Feature selection and model learning are treated as two separate steps. These procedures rely on statistical characteristics of the data such as correlation, distance and information, to select the most important features. Embedded procedures incorporate feature selection as part of the learning model, as seen in neural nets. We focus on the model-independent filter procedures for feature selection, because of their classifier independence, simplicity and computational efficiency [3]. Specifically, we consider Mutual Information (MI) based criteria for feature selection. MI is a probabilistic measure that captures the ‘correlation’ between random variables (see Figure (1)). Whereas standard correlation captures linear relationships between variables, MI can capture non-linear dependencies between variables [4].

Since our aim is to develop better classifier models using feature selection, we select the best subset of features that together have the highest MI with the class variable. Estimating MI between a subset of features requires the estimation of high dimensional joint probability distributions, which in turn necessitates exponentially large amounts of data. We circumvent this hurdle by invoking the assumption of conditional independence between the features. This reduces the MI estimation problem to a Conditional Mutual Information (CMI) estimation problem, involving three features at a time. However, selecting the subset of features with the highest CMI turns out to be a very hard combinatorial problem. We model feature selection as a Binary Quadratic Problem (BQP), which is NP-hard. We introduce approximations to solve the BQP taking inspiration from solutions to other related NP-hard problems like $k$-Sparse-PCA and Densest-$k$-Subgraph. We also evaluate our procedure by comparing it with less optimal, but computationally more efficient iterative procedures for feature selection. We evaluate our feature selection using classification accuracies across multiple datasets.

Motivation: Theoretical underpinnings of using CMI for feature selection have not been clearly outlined in literature. Also, current global feature selection methods like Semidefinite Programming (SDP) do not scale well. We therefore aim to provide better insights into the approximations and assumptions involved in CMI based feature selection, and introduce efficient approximations to solving the NP-hard problem of feature selection from related problems.

Contributions: We have demonstrated that the class posterior distribution can be approximated by selecting a subset of variables with the highest MI with the class variable. We have proved that the MI between a subset of variables and the class random variable, reduces to CMI between pairs of variables under the assumptions of conditional and class-conditional independence. We have modeled the NP-hard problem of feature selection as a Binary Quadratic Problem (BQP) and demonstrated our feature selection method across multiple datasets.

The rest of our paper is organized as follows. In Sec. II, we outline the problem of feature selection and formulate the BQP and its approximate solutions. Sec. III compares our work with existing feature selection procedures. We conclude with experiments in Sec. IV followed by discussion in Sec. V.

II Feature Selection using Conditional Mutual Information

We consider the standard classification setting where we are presented with i.i.d data $\mathscr{D}=\{(\mathbf{x}^{i},y^{i});i=1\ldots m\}$. Each data point $\mathbf{x}^{i}\in\mathbb{R}^{n}$, is regarded as an instance of a set of $n$ continuous random variables $X=\{X_{1},X_{2},\ldots,X_{n}\}$. We interchangeably use the terms features or variables to denote $X$ or its subsets. The dependent class label $y^{i}$, is considered to be an instance of a discrete random variable $Y$, that takes values in $[1,\ldots,c]$. We define $\mathbb{S}$ to be a subset of $k$ feature indices, where $1\leq k\leq n$ and $X_{\mathbb{S}}$, the subset of $k$ features indexed by $\mathbb{S}$. Likewise, $\mathbb{\tilde{S}}$ is the subset of left over indices, i.e. $\{1,\ldots,n\}\backslash\mathbb{S}$ and $X_{\mathbb{\tilde{S}}}$ is the subset of leftover $(n-k)$ features indexed by $\mathbb{\tilde{S}}$, i.e. $X\backslash X_{\mathbb{S}}$. We consider $p(Y|X)$, to be the posterior probability on the dataset $\mathscr{D}$. We now state our problem.

Problem II.1

To estimate the optimal subset $\mathbb{S}^{*}\subset\{1,\ldots,n\}$, of feature indices, such that $p(Y|X_{\mathbb{S}^{*}})$ is a good approximation to $p(Y|X)$.

Theorem II.1

$p(Y|X)=p(Y|X_{\mathbb{S}})$, if and only if Mutual Information $I(X;Y)=I(X_{\mathbb{S}};Y)$

Proof:

Joint probability $p(X,Y)$ can be factored as

[TABLE]

When viewed as elements in a Bayesian Network, the random variables can be thought to form a Markov Chain $X_{\mathbb{\tilde{S}}}\rightarrow X_{\mathbb{S}}\rightarrow Y$ where $Y$ is conditionally independent of $X_{\mathbb{\tilde{S}}}$ given $X_{\mathbb{S}}$. Therefore, the CMI $I(X_{\mathbb{\tilde{S}}};Y|X_{\mathbb{S}})=0$. With $I(X;Y)=I(X_{\mathbb{S}};Y)+I(X_{\mathbb{\tilde{S}}};Y|X_{\mathbb{S}})$, the necessary condition is proved. The statements in the proof are also true in the reverse order. ∎

Theorem(II.1) provides the justification for choosing a MI based approach to determining the optimal $\mathbb{S}^{*}$. $I(X;Y)$ depends only on $\mathscr{D}$ and is a constant that is shared between $I(X_{\mathbb{S}};Y)$ and $I(X_{\mathbb{\tilde{S}}};Y|X_{\mathbb{S}})$. It is very likely that $I(X_{\mathbb{\tilde{S}}};Y|X_{\mathbb{S}})$ is a non-zero value for every $\mathbb{S}$ except for the most trivial case with $X_{\mathbb{\tilde{S}}}=\emptyset$. For a fixed value of $k$, we therefore try to estimate the optimal subset $\mathbb{S}$ that maximizes $I(X_{\mathbb{S}};Y)$.

II-A The Binary Quadratic Problem

Given a subset of features $X_{\mathbb{S}}$, in order to evaluate $I(X_{\mathbb{S}};Y)$, we need to estimate the joint probability distribution $p(X_{\mathbb{S}},Y)$, which is often intractable. To simplify the estimation of the joint distribution, we introduce the assumption of conditional independence between features in the spirit of Naïve Bayes and [5]. We first introduce a new term $\mathbb{S}_{i}$ as the subset $\mathbb{S}$ without the index $i$ also denoted as $\mathbb{S}\backslash\{i\}$. We now outline the assumption that will simplify the joint distribution.

Assumption II.1

For a set of selected features $\{X_{\mathbb{S}_{i}}\cup X_{i}\}$, the features $X_{\mathbb{S}_{i}}$ are conditionally independent and class-conditionally independent given $X_{i}$, i.e. $p(X_{\mathbb{S}_{i}}|X_{i})=\prod_{j\in\mathbb{S}_{i}}p(X_{j}|X_{i})$ and $p(X_{\mathbb{S}_{i}}|X_{i},Y)=\prod_{j\in\mathbb{S}_{i}}p(X_{j}|X_{i},Y)$.

Theorem II.2

If Assumption II.1 is true, and $X_{i}\in X_{\mathbb{S}}$ then, $I(X_{\mathbb{S}};Y)=I(X_{i};Y)+\sum_{j\in\mathbb{S}_{i}}I(X_{j};Y|X_{i})$

Proof:

The MI between $X_{\mathbb{S}}$ and $Y$ is given by:

[TABLE]

In the above derivation, we have first applied the MI Chain Rule $I(A,B;C)=I(A;C)+I(B;C|A)$. We have then applied the MI rule $I(A;B|C)-I(A;B)=I(A;C|B)-I(A;C)$ and expressed the two trailing MI terms in terms of entropy. Applying Assumption II.1 to (2), the high order entropy terms can be replaced with summations and reduced further using $H(X_{j}|X_{i})-H(X_{j}|X_{i},Y)=I(X_{j},Y|X_{i})$, to yield,

[TABLE]

∎

Since we would like to select $\mathbb{S}$ that maximizes (3), we can formulate the global feature selection problem as,

[TABLE]

This is equivalent to the constrained Binary Quadratic problem,

[TABLE]

where, $\mathbf{Q}$ is a $[n\times n]$ non-negative matrix with $Q_{ii}$ = $I(X_{i};C)$ and $Q_{ij}$ = $I(X_{j};Y|X_{i})$ and the non-zero indices of the solution $\mathbf{x}$, constitute $\mathbb{S}$.

II-B Solving the Binary Quadratic Problem

We aim to solve the BQP, where $\mathbf{Q}$ is a symmetric and possibly indefinite matrix with non-negtiave elements, i.e., $Q_{ij}\geq 0$ for all $1\leq i,j\leq n$. Although this problem is well defined, it is highly nonconvex due to the nonconvex constraint. And it is also known that this binary quadratic problem is NP-hard [6]. So, it is difficult to find the optimal value in practice. Therefore, we aim to find an approximate solution to the BQP. Approximation methods such as linear relaxation, spectral relaxation, semidefinite programming, truncated power method and low rank bilinear approximation have been applied to solve this family of NP-hard problems. We will now briefly review some of these relaxation methods.

II-B1 Linear Relaxation

By introducing a new variable $w_{ij}$, we are able to linearize the quadratic term. With ${\mathbf{x}}=[x_{1},\ldots,x_{n}]^{\top}$, we formulate the following:

[TABLE]

LP1 is equivalent to BQP, which is also NP-hard. LP1 is simplified by relaxing $w_{ij}\in[0,1]$. One of the optimality conditions is then given by, $2w_{ij}=x_{i}+x_{j}$. This relaxation reduces LP1 to LP2 which is given by,

[TABLE]

Since $Q_{ij}\geq 0$, the maximum value for Linear is equivalent to the $k$ largest column (or row) sum of $\mathbf{Q}$. The solution to Linear guarantees a tight lower bound to the BQP (see appendix). In our work, we use the solution from Linear to initialize the input to other algorithms.

II-B2 Truncated Power Method

Truncated power (TPower) method aims to find the largest $k$-sparse eigenvector. Given a positive semidefinite matrix $A$, the largest $k$-sparse eigenvalue can be defined as follows [7]:

[TABLE]

Matrix $A$ is required to be positive semidefinte, but TPower method can be extended to deal with general symmetric matrices by setting $A\leftarrow(A+\tilde{\lambda}I_{n\times n})$ where $\tilde{\lambda}>0$ such that $(A+\tilde{\lambda}I_{n\times n})\in\mathbb{S}_{+}^{n}$. The truncated power method is given as follows. Starting from an initial $k$-sparse vector ${\mathbf{x}}_{0}$, at each iteration $t$, we multiply the vector ${\mathbf{x}}_{t-1}$ by $A$ and then truncate the entries of $A{\mathbf{x}}_{t-1}$ to zeros and set the largest $k$ entries to 1. TPower can benefit from a good starting point. We use the solution from Linear as the initial sparse vector ${\mathbf{x}}_{0}$.

II-B3 Low Rank Bilinear Approximation

The low rank bilinear approximation procedure has been applied to solve the $k$-Sparse-PCA [8] and the Densest-$k$-Subgraph [9]. It approximates the solution to a BQP by applying a bilinear relaxation which is given as,

[TABLE]

$\text{BQP}_{b}$ provides a good (2$\rho$-approximation [9]) approximation to BQP and can be solved in polynomial time using a $d$-rank approximation of ${\mathbf{Q}}$. The authors in [9] have developed the Spannogram algorithm to estimate a candidate set (termed the Spannogram $\mathscr{S}$) of vector pairs $({\mathbf{x}},{\mathbf{y}})$ with $k$ features. One of the vectors from the pair that maximizes $\text{BQP}_{b}$, is the bilinear approximate solution to BQP.

II-B4 Related BQP Approximation Methods

For the sake of completeness we briefly mention two other techniques that have been applied to approximate the BQP in the domain of feature selection. The authors in [10] propose two relaxation techniques, (1) Spectral relaxation and (2) Semidefinite Programming relaxation. In Spectral Relaxation, the constraint on the values of ${\mathbf{x}}$ are relaxed to continuous values. The values of ${\mathbf{x}}$ being positive, can be interpreted as feature weights. The solution to Spectral has been shown to be the largest eigenvector of $\mathbf{Q}$ [10]. Semidefinite Programming relaxation (SDP) is a closer approximation to the BQP than Spectral. The BQP is approximated by a trace maximization problem using semidefinite relaxation. The approximate solution ${\mathbf{x}}$, to the BQP, is obtained from the SDP solution through random projection rounding based on Cholesky decomposition [11]. For our experiments, we apply random rounding with 100 projections. For further details, please refer to [10].

III Existing Methods For MI Based Feature Selection

In this section we will discuss existing filter based Mutual Information measures for feature selection. We categorize these procedures as Greedy (iterative) Selectors and Global Selectors. We limit our discussion to the best feature selectors we found in our studies. For a broader perspective on feature selection procedures, we suggest the works of Tang et al. [12], Guyon and Elisseeff [3] and Duch [13].

III-A Greedy Feature Selectors

Greedy selectors usually begin with an empty set $\mathbb{S}$, and iteratively add to it the most important feature indices, until a fixed number of feature indices are selected or a stopping criterion is reached. MI between the random variables (features and label) provides the ranking for the features. The most basic form of the scoring function is Maximum Relevance (MaxRel) [14], where the score is simply the MI between the feature and the class variable. To account for the redundancy $I(X_{i};X_{j})$ between features, Peng et al. [15], introduced the Maximum Relevance Minimum Redundancy (MRMR) criterion, which selects features with maximum relevance to the label and minimum redundancy between each other. A greedy procedure very closely related to our technique is the Joint Mutual Information (JMI), that was developed by Yang and Moody [16], and later by Meyer et al. [17]. In our experiments (Sec. IV), we evaluate the features selected by these iterative procedures across multiple datasets.

III-B Global Feature Selectors

There is limited work on MI based global feature selection. In [4], Rodriguez-Lujan et al. introduced the Quadratic Programming Feature Selection (QPFS). This method can be viewed as a global alternative to MRMR. The second global technique proposed by Nguyen et al. [10], is related to our method presented in (BQP). Nguyen et al. [10], model a global feature selection problem based on the CMI matrix ${\mathbf{Q}}$, and apply Spectral and SDP methods to approximate the solution.

IV Experiments

We consider two factors to evaluate the feature selection algorithms discussed so far, viz., time complexity and classification accuracies. For the global methods, since we are approximating the solution, we also consider the tightness of the approximation. We conducted our experiments using MATLAB on an Intel Core-i7 2.3 GHz processor with 16GB of memory.

IV-A Feature Selectors: a test of scalability

Greedy feature selectors have time complexities of the order $O(nk)$, which is negligible compared to the time complexities of the global feature selectors. Table I lists the time complexities for the global algorithms. To study time complexities, we conduct multiple experiments $(n,k)$, where we simulate a CMI matrix ${\mathbf{Q}}$ by a random positive symmetric matrix of size $[n\times n]$, and select $k$ features. The time complexity for experiment $(n,k)$, is the average time of convergence over 10 runs. We use the same set of random matrices for each of the algorithms in the experiment. Figure (2) depicts the convergence times for different experiments. Linear algorithm is the most efficient, followed closely by Spectral and TPower methods. We used the CVX [18] implementation with SDPT3 solver [19], for all our SDP experiments. The SDP solver has a huge memory footprint and with matrix sizes $n\geq 700$, we run into ‘Out of Memory’ errors. For the LowRank method, we used the following parameters, $d=3,\epsilon=0.1,\delta=0.1$ for all our experiments. Please refer to [9] for more details on the LowRank.

IV-B BQP Methods: a test of approximation

For the next set of experiments, we evaluated the degree of approximation for the global algorithms. Since we do not know the optimal solution for the BQP, we compared the methods on their relative objective values. We estimated the binary feature vector $\mathbf{x}$ after applying each of the methods and then evaluated the objective value $\mathbf{x}^{\top}\mathbf{Q}\mathbf{x}$. We evaluated the percentage difference of every algorithm’s objective value with the Linear method’s objective. Similar to the experimental evaluation used for time complexity, we generated random data for each experiment $(n,k)$ and averaged the values over 10 random runs. Figure (3) presents the results of the experiment. TPower and LowRank displayed the largest percentage increase from the Linear. Since TPower and LowRank approximate the BQP better than other methods, they must therefore be better feature selectors compared to Linear, Spectral and SDP. The TPower is also very efficient in terms of execution time and would be the ideal feature selector when considering both speed and accuracy.

IV-C Feature Selectors: a test of classification error

In this section we compare the TPower and the LowRank with other algorithms in terms of classification accuracies. For our experiments, we chose 13 publicly available datasets that are widely used to study MI based feature selection as in [4, 5, 10, 15]. The details of the datasets are captured in Table II. We performed feature selection for a set of $k$ values and estimated the classifier performance across all values of $k$. Starting at $k=10$, we incremented in steps of $1$ till $n$ or $100$, whichever is smaller. We evaluated the classifier performance using Leave-One-Out cross validation (if $m\leq 100$) or 10-fold cross validation and obtained the cross validation errors(%) for each fold. Since the average error across all values of $k$ is not a good measure of the classifier performance, we applied the paired t-test, as also mentioned in [4, 10, 22], across the cross validation folds. For a fixed dataset and a fixed value of $k$, to compare TPower with say, MaxRel, we applied the one sided paired t-test at 5% significance over the error(%) of the cross validation folds for the two algorithms. We set the performance of TPower vs. MaxRel to win = $w$, tie = $t$ and loss = $l$, based on the largest number of t-test decisions across all the $k$ values. Along the lines of earlier studies in feature selection, we used linear SVM as the classifier. To estimate the CMI we need to discretize the features. We believe the role of discretization is not unduly critical as long as it is consistent across all the experiments. We discretized the features using the Class Attribute Interdependence Maximization (CAIM) algorithm developed by Kurgan and Cios et al. [23]. Feature selection was performed on discretized data but the classification (after feature selection) was performed on the original feature space. Using the above procedure, we compared the performance of TPower and LowRank with all the other algorithms. Tables III and IV display the results of the experiment. The values in Table III (likewise IV) correspond to the difference in the average of classification error(%) between TPower (likewise LowRank) and all the other algorithms. From the results in these tables, we find that TPower and LowRank do well on most of the datasets across all the algorithms. When compared against each other LowRank does better than TPower. For the sake of brevity we have not displayed the comparisons between other pairs of algorithms. The win/tie/loss numbers by themselves do not provide a complete picture of the comparison. The difference in the average error also needs to be taken into account to assess the performance. A large percentage of negative values in the columns and their magnitudes indicate the low error values in classification for TPower and LowRank. Figure (4) displays the average classification error(%) trends for varying values of $k$ for $3$ datasets. Figures (4(a)) and (4(d)) for the Colon dataset, suggest that the addition of more features does not necessarily reduce classification error. Classification error trends also help us cross validate the best value of $k$ for a dataset. For a given dataset, the error trends between the global and greedy procedures follow a similar pattern. This perhaps indicates that nearly similar features are being selected using both types of methods. We also note that for huge datasets with large values of $n$, greedy methods may not be a bad choice for feature selection.

V Discussion and Conclusions

Feature selection is a NP-hard problem and newer methods to approximate the solution will help drive research in this area. We have demonstrated that current methods applying MI and CMI for feature selection, inherently assume conditional independence between features. The conditional independence assumption limits the number of features compared to only 2 or 3 features at a time. There is need to derive better measures to approximate the importance of a group of selected features. To estimate the probability distributions, we had to discretize the features. We have not studied the effect of discretization in our work. Progress along all of these fronts will provide directions to improve MI based feature selection. In conclusion, we can state that both TPower and LowRank perform better than existing global and iterative techniques across most of the datasets. While LowRank slightly outperforms TPower, it does not compare well with regards to time.

VI Acknowledgments

This material is based upon work supported by the National Science Foundation (NSF) under Grant No:1116360. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF.

We express our thanks to Nguyen X. Vinh [10], for providing his implementation for comparison.

Appendix A LP2 provides a Good Lower bound for BQP

We would like to study the goodness of approximation to BQP provided by the solution to LP2. We define an equivalent problem in terms of BQP.

Proposition A.1

[TABLE]

is equilvalent to BQP.

Proof:

$||{\mathbf{Q}}||_{1}$ is the sum of all elements in ${\mathbf{Q}}$. Since ${\mathbf{Q}}_{ij}\geq 0$, BQP $\leq||{\mathbf{Q}}||_{1}~{}\forall{\mathbf{x}}$. Therefore, $q({\mathbf{x}})\leq 0,~{}\forall{\mathbf{x}}$. Since $||{\mathbf{Q}}||_{1}$ is a constant for a matrix, under the same set of constraints,

[TABLE]

∎

We define some new quantities for the derivation of the bound. Let ${\mathbf{x}}^{*}$ be the solution of BQP. Let $\bar{{\mathbf{x}}}$ be the solution of LP2. Since $||{\mathbf{Q}}||_{1}=\sum_{ij}{\mathbf{Q}}_{ij}$, we can expand $||{\mathbf{Q}}||_{1}$ in terms of any binary vector ${\mathbf{x}}$. Specifically we define $||{\mathbf{Q}}||_{1},$ in terms of $\bar{{\mathbf{x}}}$,

Definition A.1

[TABLE]

Lemma A.1

$||{\mathbf{Q}}||_{1}-\bar{{\mathbf{x}}}^{\top}{\mathbf{Q}}\bar{{\mathbf{x}}}\geq Q^{1}$**

Proof:

From (7) we have,

[TABLE]

∎

Let $Q^{*}$ denote the maximum value of BQP and let $Q^{*}_{LP2}$ denote maximum value of LP2. If ${\mathbf{x}}^{*}$ is the solution of BQP and $\bar{{\mathbf{x}}}$ is the solution of LP2. We have the following result:

Lemma A.2

$Q^{*}_{LP2}\geq Q^{*}$**

Proof:

[TABLE]

∎

We are now ready to state the bound for LP2.

Theorem A.1

[TABLE]

Proof:

From Lemma (A.2), we have:

[TABLE]

∎

The last statement (16) is arrived at by adding $-2||{\mathbf{Q}}||_{1}$ on both sides. Since $q({\mathbf{x}})\leq 0$, $2q({\mathbf{x}}^{*})\leq q(\bar{{\mathbf{x}}})$ implies that $q(\bar{{\mathbf{x}}})$ is a lower bound for $q({\mathbf{x}}^{*})$. We are therefore guaranteed a lower bound for QBP by solving LP2 and (16) provides the tightness of the bound.