High Dimensional Classification through $\ell_0$-Penalized Empirical Risk Minimization

High Dimensional Classification through $_0$-Penalized Empirical Risk Minimization

Introduction

Binary classification is concerned with learning a binary classifier that can be used to categorize objects into one of two predefined statuses. It arises in a wide range of applications and has been extensively studied in the statistics and machine learning literature. For comprehensive surveys and discussions on the binary classification methods, see e.g. DGL:1996, vapnik2000, Lugosi:2002, boucheron2005 and hastie2009. Solving for the optimal binary classifier by minimizing the empirical misclassification risk is known as an empirical risk minimization (ERM) problem. There has been a massive research interest in high dimensional classification problems where the dimension of the feature vector used to classify the object's label can be comparable with or even larger than the available training sample size. It is known (see e.g., bickel2004 and fan2008) that working directly with a high dimensional feature space can result in poor classification performance. To overcome this, it is often assumed that only a small subset of features are important for classification and feature selection is performed to mitigate the high dimensionality problem. See fan2011 for an overview on the issues and methods for high dimensional classification.

In this paper, we study the ERM based binary classification in the setting with high dimensional vectors of features. We propose an $\ell _{0}$-penalized ERM procedure for classification by minimizing over a class of linear classifiers the empirical misclassification risk with a penalty on the number of selected features. Here, the $\ell _{0}$-norm of a real vector refers to the number of non-zero components of the vector. When the Bayes classifier, which is the optimal classifier that minimizes the population misclassification risk, is also of the linear classifier form and respects a sparsity condition, we show that this penalized ERM classification approach can yield a sparse solution for feature selection with high probability. Moreover, we derive non-asymptotic bound on the excess misclassification risk and establish its rate of convergence.

There are alternative ERM based approaches to the high dimensional binary classification problem. Greenshtein2006, Jiang:10 and chen2018 studied the best subset variable selection approach where the ERM problem is solved subject to a constraint on a pre-specified maximal number of selected features. Jiang:10 further showed that the best subset ERM problem can be approximated by the $\ell _{0}$-penalized ERM problem. However, they did not establish theoretical results characterizing the size of the subset of feature variables selected under the $\ell _{0}$-penalized ERM approach. Neither did they provide numerical algorithms for solving the $\ell _{0}$-penalized estimation problem.

In the present paper, we take the $\ell _{0}$-penalized ERM approach and develop a computational method for implementation. We show that there is a high probability in large samples that the resulting number of selected features under this penalized estimation approach can be capped above by an upper bound which can be made arbitrarily close to the unknown smallest number of features that are relevant for classification. Our penalized ERM approach is also closely related to the method of structural risk minimization (see, e.g., DGL:1996) where the best classifier is selected by solving a sequence of penalized ERM problems over an increasing sequence of spaces of classifiers with the penalty depending on the complexity of the classifier space measured in terms of the Vapnik-Chervonenkis (VC) dimension. As will be discussed later, our approach can also be interpreted in a similar fashion yet with a different type of complexity penalty.

For implementation, we show that the $\ell _{0}$-penalized ERM problem of this paper can be equivalently reformulated as a mixed integer linear programming (MILP) problem. This reformulation enables us to employ modern efficient MIO solvers to solve our penalized ERM problem. Well-known numerical solvers such as CPLEX and Gurobi can be used to effectively solve large-scale MILP problems. See Nemhauser1999 and bertsimas2005 for classic texts on the MIO theory and applications. See also junger2009, achterberg2013 and bertsimas2016 for discussions on computational advances in solving the MIO problems.

The present paper is organized as follows. In Section (ref), we describe the binary classification problem and set forth the $\ell _{0}$-penalized ERM approach. In Section (ref), we establish theoretical properties of the proposed classification approach. In Section (ref), we provide a computational method using the MIO approach. In Section (ref), we conduct a simulation study on the performance of the $\ell _{0}$-penalized ERM approach in high dimensional binary classification problems. We then conclude the paper in Section (ref). Proofs of all theoretical results of the paper are collated in Appendix (ref).

An $\ell_0$-Penalized ERM Approach

Let $Y\in \{0,1\}$ be the binary label or outcome of an object and $X$ a $ \left( p+1\right) $ dimensional feature vector of that object. Write $ X=(X_{1},\widetilde{X})$, where $X_{1}$ is a scalar random variable that is always included and has a positive effect and $\widetilde{X}$ is the $p$ dimensional subvector of $X$ subject to feature selection. For $x\in \mathcal{X}$, let

where $\mathcal{X}$ is the support of $X$, $\theta $ is a vector of parameters, and $1\left\{ \cdot \right\} $ is an indicator function that takes value 1 if its argument is true and 0 otherwise.

We consider binary classification using linear classifiers of the form ((ref)). Since the condition $X_{1}+\widetilde{X}^{\prime }\theta \geq 0 $ is invariant with respect to any positive scalar that multiplies both sides of this inequality, working with the classifier ((ref)) amounts to normalizing the scale by setting the coefficient of $X_{1}$ to be unity. For any $p$ dimensional real vector $\theta $, let $\left\Vert \theta \right\Vert _{0}\equiv \sum\nolimits_{j=1}^{p}1\{\theta _{j}\neq 0\}$ be the $\ell _{0}$-norm of $\theta $. Assume that the researcher has a training sample of $n$ independent identically distributed (i.i.d.) observations $ \left( Y_{i},X_{i}\right) _{i=1}^{n}$ of $(Y,X)$. We allow the dimension $p$ to be potentially much larger than the sample size $n$. We estimate the coefficient vector $\theta $ by solving the following $\ell _{0}$-penalized minimization problem

where $\Theta \subset \mathbb{R}^{p}$ denotes the parameter space, and, for any indicator function $b:\mathcal{X}\mapsto \left\{ 0,1\right\} $,

and $\lambda $ is a given non-negative tuning parameter of the penalized minimization problem.

The function $S_{n}(b)$ is known as the empirical misclassification risk for the binary classifier $b$. Minimization of $S_{n}(b)$ over the class of binary classifiers given by ((ref)) is known as an empirical risk minimization (ERM) problem. The penalized ERM approach ((ref)) enforces dimension reduction by attaching a higher penalty to a classifier $ b_{\theta }$ which uses more object features for classification. Let $ \widehat{\theta }$ be a solution to the minimization problem ((ref)). We shall refer to the resulting classifier ((ref)) evaluated at $\widehat{\theta }$ as an $\ell _{0}$-penalized ERM classifier.

For any $m\geq 0$, let

where

That is, $\mathcal{B}_{m}$ is the class of all linear classifiers in ((ref)) whose $\theta $ vector has no more than $m$ non-zero components. For $m\in \{0,1,...,p\}$, let

Then it is straightforward to see that the minimized objective value of the penalized ERM problem ((ref)) is equivalent to that of the problem

In other words, our approach is akin to the method of structural risk minimization as it amounts to solving ERM problems over an increasing sequence of classifier spaces $\mathcal{B}_{m}$ which carries a complexity penalty $\lambda m$. In the next section, we will set forth regularity conditions on the penalty tuning parameter $\lambda $ and establish theoretical properties for our classification approach.

Theoretical Properties

In this section, we study theoretical properties of the $\ell _{0}$-penalized ERM classification approach. Let $F$ denote the joint distribution of $\left( Y,X\right) $. For any indicator function $b:\mathcal{X}\mapsto \left\{ 0,1\right\} $, let

For $x\in \mathcal{X}$, let

For any measurable function $f:\mathcal{W\mapsto }\mathbb{R}$, let $ \left\Vert f\right\Vert _{1}=E\left[ \left\vert f(W)\right\vert \right] $ denote the $L_{1}$-norm of $f$. The functions $\eta $ and $b^{\ast }$ as well as the $L_{1}$-norm $\left\Vert \cdot \right\Vert _{1}$ depend on the data generating distribution $F$. It is straightforward to see that, for any binary classifier $b$,

so that $S(b)$ is minimized at $b=b^{\ast }$. The optimal classifier $ b^{\ast }$ is known as the Bayes classifier in the classification literature.

We assess the predictive performance of the $\ell _{0}$-penalized ERM approach by bounding the excess risk

The difference $U_{n}$ is non-negative by ((ref)). Hence, a good classifier will result in a small value of $U_{n}$ with a high probability and also on average.

We impose the following assumption.

Let $q$ denote the smallest value of non-negative integers $d$ satisfying $ b^{\ast }\in \mathcal{B}_{d}$. By Condition (ref), such $q$ value is finite and always exists. Condition (ref) implies that the Bayes classifier $b^{\ast }$ admits a linear threshold crossing structure in the sense that the equivalence

holds almost surely for some $\theta \in \Theta _{q}$, where $q$ can be interpreted as the sparsity parameter associated with $b^{\ast }$, which is unknown in this binary classification problem. Moreover, the assumption that $q\leq p$ implies that the feature vector $(X_1, \widetilde{X})$ is rich enough to embody those relevant ones for constructing the Bayes classifier.

For any two real numbers $x$ and $y$, let $x\vee y\equiv \max \{x,y\}$ and $ x\wedge y\equiv \min \{x,y\}$. For any $x\geq 0$, let $\left\lceil x\right\rceil $ and $\left\lfloor x\right\rfloor $ respectively denote the integer ceiling and floor of $x$. We impose the following condition on the growing rates of $\lambda $ relative to the sample size.

Let

The estimate $\left\Vert \widehat{\theta }\right\Vert _{0}$ corresponds to the number of features selected under the $\ell _{0}$-penalized ERM approach. We now provide a result on the statistical behavior of $\left\Vert \widehat{\theta }\right\Vert _{0}$, which sheds lights on the dimension reduction performance of our penalized estimation approach.

For any fixed $\epsilon \in (0,1)$, we can deduce from Theorem (ref) that $P\left( \left\Vert \widehat{\theta }\right\Vert _{0}>s\right) \longrightarrow 0$ as $p\vee n\longrightarrow \infty $. Moreover, this theorem implies that our approach is effective in reducing the feature dimension in the sense that, with high probability in large samples, the number of selected features is capped above by the quantity ((ref)), which can be made arbitrarily close to the true sparsity $q$ in the classification problem. Specifically, if $\epsilon $ turns out to be smaller than $1/\left( q+1\right) $, the result ((ref)) implies that $P\left( \left\Vert \widehat{\theta }\right\Vert _{0}>q+1\right) $ tends to zero exponentially in $r_{n}$.

The next theorem characterizes the predictive performance of the $\ell _{0}$ -penalized ERM approach.

Theorem (ref) implies that the tail probability of $U_{n}$ decays to zero exponentially in $r_{n}$. Moreover, inequality ((ref)) together with the fact that $U_{n}\leq 1$ immediately implies that

By Condition (ref) and ((ref)), we can therefore deduce that

which converges to zero whenever

The rate condition ((ref)) allows the case that

In other words, the $\ell _{0}$-penalized ERM classification approach is risk-consistent even when the dimension of the input feature space ($p$) grows exponentially in sample size, provided that the number of truly effective features ($q$) can only grow at a polynomial rate.

We shall provide some further remarks on the convergence rate result ((ref)). Condition (ref) implies that the space $ \mathcal{B}_{q}$ contains the Bayes classifier $b^{\ast }$. Thus, if the value of $q$ were known, one could performed classification via the $\ell _{0}$-constrained ERM approach where the empirical risk $S_{n}(b)$ is minimized with respect to $b\in \mathcal{B}_{q}$. The lower bound on the VC dimension of the classifier space $\mathcal{B}_{q}$ grows at rate $O(q\ln \left( p/q\right) )$ (abramovich2017). Hence, the rate $O( \sqrt{n^{-1}q\ln \left( p/q\right) })$ is the minimax optimal rate at which the excess risk converges to zero under this constrained estimation approach (DGL:1996). In view of this, suppose $p$ grows at a polynomial or exponential rate in $n$. Then our rate result ((ref)) is nearly oracle in the sense that, when $q$ grows at rate $O(\ln n)$, the rate ((ref)) remains close within some $\ln n$ factor to the optimal rate attained under the case of known $q$. Moreover, both rates coincide and reduce to $O(\sqrt{n^{-1}\ln p} ) $ when the value of $q$ does not increase with the sample size.

Computational Algorithms

While the ERM approach to binary classification is theoretically sound, its implementation is computationally challenging and is known to be an NP (Non-deterministic Polynomial time) hard problem (Johnson:78). Florios:Skouras:08 developed a mixed integer optimization (MIO) based computational method and provided numerical evidence demonstrating effectiveness of the MIO approach to solving the ERM type optimization problems. kitagawa2018 and mbakop2018 adopted the MIO solution approach to solving the optimal treatment assignment problem which is closely related to the ERM based classification problem. The MIO approach is also useful for solving problems of variable selection through $\ell _{0}$ -norm constraints. See bertsimas2016 and chen2018 who proposed MIO based computational algorithms to solving the $\ell _{0}$-constrained regression and classification problems respectively.

Motivated by these previous works, we now present an MIO based computational method for solving the $\ell _{0}$-penalized ERM problem. Given that $ Y_{i}\in \{0,1\}$, solving the problem ((ref)) amounts to solving

We assume that the parameter space $\Theta $ takes the form

where $\underline{c}_{j}$ and $\overline{c}_{j}$ are lower and upper parameter bounds such that $-\infty <\underline{c}_{j}\leq \theta _{j}\leq \overline{c}_{j}<\infty $ for $j\in \{1,...,p\}$.

Our implementation builds on the method of mixed integer optimization. Specifically, we note that the minimization problem ((ref)) can be equivalently reformulated as the following mixed integer linear programming problem:

where $\delta $ is a given small and positive real scalar (e.g. $\delta =10^{-6}$ as in our numerical study), and

We now explain the equivalence between ((ref)) and ((ref)). Given $\theta $, the inequality constraints ((ref)) and the dichotomization constraints ((ref)) enforce that $ d_{i}=1\{X_{1i}+\widetilde{X}_{i}^{\prime }\theta \geq 0\}$ for $i\in \{1,...n\}$. Moreover, the on-off constraints ((ref)) and ((ref)) ensure that, whenever $e_{j}=0$, the value $\theta _{j}$ must also be zero so that the $j$th component of the feature vector $\widetilde{X}$ is excluded in the resulting $\ell _{0}$-penalized ERM classifier. The sum $\sum\nolimits_{j=1}^{p}e_{j}$ thus captures the number of non-zero components of the vector $\theta $. As a result, both minimization problems ((ref)) and ((ref)) are equivalent. This equivalence enables us to employ modern MIO solvers to solve for $\ell _{0}$-penalized ERM classifiers. For implementation, note that the values $\left( M_{i}\right) _{i=1}^{n}$ in the inequality constraints ((ref)) can be computed by formulating the maximization problem in ((ref)) as linear programming problems, which can be efficiently solved by modern numerical solvers. Hence these values can be easily computed and stored as inputs to the MILP problem ((ref)).

Simulation Study

In this section, we conduct simulation experiments to study the performance of our approach. We consider a simulation setup similar to that of chen2018 and use the following data generating design. Let $V=(V_{1},...,V_{p})$ be a multivariate normal random vector with mean zero and covariance matrix $ \Sigma $ with its element $\Sigma _{i,j}=\left( 0.25\right) ^{\left\vert i-j\right\vert }$. The binary outcome is generated according to the following specification:

where $\theta ^{\ast }$ denotes the true data generating parameter value, $ X=(X_{1},\widetilde{X})$ is a $(p+1)$ dimensional feature vector with $ X_{1}=V_{1}$ and $\widetilde{X}=(1,V_{2},...,V_{p})$, and $\xi $ is a random variate that is independent of $V$ and follows the standard logistic distribution. The constant term in $\widetilde{X}$ is included to capture the regression intercept. We set $\theta _{1}^{\ast }=0$ and $\theta _{j}^{\ast }=0$ for $j\in \{3,...,p\}.$ The coefficient $\theta _{2}^{\ast }$ is chosen to be non-zero such that, among all the feature variables in $ \widetilde{X}$, only the variable $\widetilde{X}_{2}=V_{2}$ is relevant in the data generating processes (DGP). We consider the following two specifications for $\theta _{2}^{\ast }$ and $\sigma (X)$:

We used $100$ simulation repetitions in each Monte Carlo experiment. For each simulation repetition, we generated a training sample of $n=100$ observations for estimating the coefficients $\theta $ and a validation sample of $5000$ observations for evaluating the out-of-sample classification performance. We considered simulation configurations with $ p\in \{10,200\}$ to assess the classifier's performance in both the low and high dimensional binary classification problems.

We specified the parameter space $\Theta $ to be $[-10,10]^{p}$ for the MIO computation of the $\ell _{0}$-penalized ERM classifiers. Throughout this paper, we used the MATLAB implementation of the Gurobi Optimizer to solve the MIO problems ((ref)). Moreover, all numerical computations were done on a desktop PC (Windows 7) equipped with 32 GB RAM and a CPU processor (Intel i7-5930K) of 3.5 GHz. To reduce computation cost of solving the $\ell _{0}$ -penalized ERM problems, we set the MIO solver time limit to be one hour beyond which we forced the solver to stop early and used the best discovered feasible solution to construct the resulting $\ell _{0}$-penalized ERM classifier. For implementation, it remains to specify an exact form of the penalty parameter $\lambda $ in ((ref)). We set

where $v$ is a tuning constant which remains to be calibrated. The form ((ref)) implies that the value $c$ in Condition (ref) is taken to be $v\ln \ln (p\vee n)$, which will satisfy inequality ((ref)) and hence validate the probability bound ( (ref)) when $p\vee n$ is sufficiently large. Moreover, by the risk upper bound ((ref) ), the convergence rate result ((ref)) continues to hold up to a factor of $\ln \ln \left( p\vee n\right) $.

For practical applications, we recommend calibrating the tuning scalar $v$ via the method of cross validation. Yet, for this simulation study, we used a simple heuristic rule and set

The value $h$ in ((ref)) can also be computed via the MIO approach by simply removing from the MIO problem ((ref)) the constraints ((ref)) and ((ref)) as well as the binary controls $\left( e_{1},..,e_{p}\right) $ and the penalty part in the objective function. This computation is much faster as it is concerned with one-dimensional optimization. The rationale behind the choice ((ref)) is as follows. Note that ((ref)) corresponds to an ERM classification using classifier ((ref)) where $\widetilde{X}$ only consists of the intercept term, and ((ref)) corresponds to an estimate of the variance of the misclassification loss under such a simple classification rule. Intuitively speaking, the value $v$ captures the variability of the empirical risk under a parsimonious feature space specification. From the bias and variance tradeoff perspective, when this variability is small, we may as well increase the classifier flexibility by attaching a small penalty in the penalized ERM procedure so as to induce a richer set of selected features for classification.

Let logit_lasso denote the $\ell _{1}$-penalized logistic regression approach friedman2010. The logit_lasso estimation approach is a computationally attractive approach that can be used to estimate high dimensional binary response models. We compared in simulations the performance of our method to that of the logit_lasso approach. We used the MATLAB implementation of the well known glmnet computational algorithms Qian2013 for solving the logit_lasso estimation problems. We did not penalize the coefficient of the feature variable $X_{1}$ so that, as in the simulations of the $\ell _{0}$-penalized ERM approach, this variable would always be included in the resulting classifier constructed under the logit_lasso approach. We calibrated the lasso penalty parameter value over a sequence of 100 values via the 10-fold cross validation procedure. We used the default setup of glmnet for constructing this tuning sequence among which we reported results based on the following two choices, $\left\{ \lambda _{opt}^{lasso},\lambda _{1se}^{lasso}\right\} $ , of the penalty parameter value. The value $\lambda _{opt}^{lasso}$ refers to the lasso penalty parameter value that minimized the cross validated misclassification risk, whereas $\lambda _{1se}^{lasso}$ denotes the largest penalty parameter value whose corresponding cross validated misclassification risk still falls within the one standard error of the cross validated misclassification risk evaluated at $\lambda _{opt}^{lasso}$ . The choice $\lambda _{1se}^{lasso}$ induces a more parsimonious estimating model and is known as the "one-standard-error" rule, which is also commonly employed in the statistical learning literature hastie2009.

We considered the following performance measures. Let $\widehat{\theta }$ denote the estimated coefficients under a given classification approach. For the logit_lasso approach, we derived $\widehat{\theta }$ by dividing the lasso-penalized logistic regression coefficients of the variables $ \widetilde{X}$ by the magnitude of that of the variable $X_{1}$. We can easily deduce that $b_{\theta ^{\ast }}(X)=1\{X_{1}+\widetilde{X}^{\prime }\theta ^{\ast }\geq 0\}$ is the Bayes classifier in this simulation design. To assess the classification performance, we report the relative risk, which is the ratio of the misclassification risk evaluated at the classifier $b_{ \widehat{\theta }}$ over that evaluated at the Bayes classifier. In each simulation repetition, we approximated the out-of-sample misclassification risk using the generated validation sample. Let $in\_RR$ and $out\_RR$ respectively denote the average of in-sample and that of out-of-sample relative risks over all the simulation repetitions.

We also examine the feature selection performance of the classification method. We say that a feature variable $\widetilde{X}_{j}$ is effectively selected if and only if the magnitude of $\widehat{\theta }_{j}$ is larger than a small tolerance level (e.g., $10^{-6}$ as used in our numerical study) which is distinct from zero in numerical computation. Let $Corr\_sel$ be the proportion of the variable $\widetilde{X}_{2}$ being effectively selected. Let $Orac\_sel$ be the proportion of obtaining an oracle feature selection outcome where the variable $\widetilde{X}_{2}$ was the only one that was effectively selected among all the variables in $\widetilde{X}$. Let $Num\_irrel$ denote the average number of effectively selected features whose true DGP coefficients are zero.

Simulation Results

We now present in Tables (ref) and (ref) the simulation results under the setups of DGP(i) and DGP(ii) respectively. From these two tables, we find that, regarding the in-sample classification performance, our method outperformed the two logit_lasso based approaches across almost all the DGP configurations in the simulation. For the out-of-sample classification performance, we see that the $\ell _{0}$ -penalized ERM classifier dominated the logit_lasso classifiers across all simulation scenarios and this dominance was more evident in the high dimensional setup with $p=200$.

Concerning the feature selection results, both Tables (ref) and (ref) indicate that all the three classification approaches had high $ Corr\_sel$ rates and hence were effective for selecting the relevant variable $\widetilde{X}_{2}$. However, the good performance in the $ Corr\_sel $ criterion might just be a consequence of overfitting, which may result in excessive selection of irrelevant variables and thus adversely impact on the out-of-sample classification performance. From the results on the $Num\_irrel $ performance measure, we note that the numbers of irrelevant variables selected under the two logit_lasso based approaches remained quite large relatively to those under the $\ell _{0}$-penalized ERM approach even though all these approaches exhibited the effect of shrinking the feature space dimension. In fact, we observe non-zero and high values of $Orac\_Sel$ for the $\ell _{0}$-classifier across all the simulation setups whereas the two logit_lasso classifiers could not induce any oracle variable selection outcome in the simulation. These feature selection performance results help to explain that the risk performance dominance of the $\ell _{0}$-penalized ERM approach could be observed even in the DGP(i) simulations where the logistic regression model was correctly specified.

Conclusions

In this paper, we study the binary classification problem in a setting with high dimensional vectors of features. We construct a binary classification procedure by minimizing the empirical misclassification risk with a penalty on the number of selected features. We establish a finite-sample probability bound showing that this classification approach can yield a sparse solution for feature selection with high probability. We also conduct non-asymptotic analysis on the excess misclassification risk and establish its rate of convergence. For implementation, we show that the penalized empirical risk minimization problem can be solved via the method of mixed integer linear programming.

There are a few topics one may consider as possible extensions. First, it might be fruitful to explore an $\ell_0$-penalized approach for regression and other estimation problems. Second, our proposed method is suitable for training samples with small or moderate size. It would be a natural step to develop a divide-and-conquer algorithm for a large-scale problem Shi-et-al:2018. Third, our approach might be applicable for developing sparse policy learning rules athey2017efficient. These are topics for further research.