Sparse group lasso and high dimensional multinomial classification.

Consider a classification problem with $K$ classes, $N$ samples, and $p$ features. Assume given a data set $(x_{1},y_{1}),\dots,(x_{N},y_{N})$ where, for all $i=1,\dots,N$, $x_{i}\in\mathbb{R}^{p}$ is the observed feature vector and $y_{i}\in\{1,\dots,K\}$ is the categorical response. We organize the feature vectors in the $N\times p$ design matrix

As in [10] we use a symmetric parametrization of the multinomial model. With $h:\{1,\dots,K\}\times\mathbb{R}^{p}\to\mathbb{R}$ given by

the multinomial model is specified by

The model parameters are organized in the $K$-dimensional vector, ${\beta}^{(0)}$, of intercept parameters together with the $K\times p$ matrix

where ${\beta}^{(i)}\in\mathbb{R}^{K}$ are the parameters associated with the $i$’th feature. The lack of identifiability in this parametrization is in practice dealt with by the penalization.

The log-likelihood is

Our interest is the sparse group lasso penalized maximum likelihood estimator. That is, $({\beta}^{(0)},\beta)$ is estimated as a minimizer of the sparse group lasso penalized negative-log-likelihood:

In our applications we let $\gamma_{J}=\sqrt{K}$ for all $J=1,\dots,p$ and $\xi_{i}=1$ for all $i=1,\dots,Kp$, but other choices are possible in the implementation. Note that the parameter grouping, as part of the penalty specification, is given in terms of the columns in (2), i.e. $m=p$.

The data sets were preprocessed before applying the multinomial sparse group lasso estimator. Two preprocessing schemes were used: normalization and standardization. Normalization entails centering and scaling of the samples in order to obtain a design matrix with row means of 0 and row variance of 1. Standardization involves centering and scaling the features to obtain a design matrix with column means of 0 and column variance of 1. Note that the order in which normalization and standardization are applied matters.

The purpose of normalization is to remove technical (non-biological) variation. A range of different normalization procedures exist for biological data. Sample centering and scaling is one of the simpler procedures. We use this simple normalization procedure for the two biological data sets in this paper. Normalization is done before and independent from the sparse group lasso algorithm.

The purpose of standardization is to create a common scale for features. This ensures that differences in scale will not influence the penalty and thus the variable selection. Standardization is an option for the sparse group lasso implementation, and it is applied as the last preprocessing step for all three example data sets.

We want to compare the performance of the multinomial sparse group lasso estimator for different values of the regularization parameter $\alpha$. Applying the multinomial sparse group lasso estimator with a given $\alpha\in[0,1]$ and $\lambda$-sequence, $\lambda_{1},\ldots,\lambda_{d}>0$, results in a sequence of estimated models with parameters $\{\hat{\beta}(\lambda_{i},\alpha)\}_{i=1,\dots,d}$. The generalization error can be estimated by cross validation [11]. For our applications we keep the sample ratio between classes in the cross validation subsets approximately fixed to that of the entire data set. Hence, we may compute a sequence, $\{\widehat{\operatorname*{Err}}(\lambda_{i},\alpha)\}_{i=1,\dots,d}$, of estimated expected generalization errors for the sequence of models. However, for given $\alpha_{1}$ and $\alpha_{2}$ we cannot simply compare $\widehat{\operatorname*{Err}}(\lambda_{i},\alpha_{1})$ and $\widehat{\operatorname*{Err}}(\lambda_{i},\alpha_{2})$, since the $\lambda_{i}$ value is scaled differently for different values of $\alpha$. We will instead compare the models with the same number of non-zero parameters and the same number of non-zero parameter groups, respectively. Define

with $\hat{\beta}(\lambda,\alpha)$ the estimator of $\beta$ for the given values of $\lambda$ and $\alpha$. That is, $\hat{\Theta}(\lambda,\alpha)$ is the number of non-zero parameter blocks in the fitted model. Note that there is a one-to-one correspondence between parameter blocks and features in the design matrix. Furthermore, we define the total number of non-zero parameters as

In particular, we want to compare the fitted models with the same number of parameter blocks. There may, however, be more than one $\lambda$-value corresponding to a given value of $\hat{\Theta}$. Thus we compare the models on a subsequence of the $\lambda$-sequence. With $\theta_{1}<\dots<\theta_{d^{\prime}}$ for $d^{\prime}\leq d$ denoting the different elements of the set $\{\hat{\Theta}(\lambda_{i},\alpha)\}_{i=1,\ldots,d}$ in increasing order we define

We then compare the characteristics of the multinomial sparse group lasso estimators for different $\alpha$ values by comparing the estimates

In this section we investigate the characteristics of the sparse group lasso estimator on simulated data sets. We are primarily interested in trends in the generalization error as $\alpha$ is varied and $\hat{\lambda}$ is selected by cross validation on a relatively small training set. We suspect that this trend will depend on the distribution of the data. We restrict our attention to multiclass data where the distribution of the features given the class is Gaussian. Loosely speaking, we suspect that if the differences in the data distributions are very sparse, i.e. the centers of the Gaussian distributions are mostly identical across classes, the lasso will produce models with the lowest generalization error. If the data distribution is sparse, but not very sparse, then the optimal $\alpha$ is in the interval $(0,1)$. For a dense distribution, typically with differences being among all classes, we expect the group lasso to perform best. The simulation study confirms this.

The mathematical formulation is as follows. Let

where $\mu_{i}\in\mathbb{R}^{p}$ for $i=1,\dots,K$ and $p=p_{a}+p_{b}$. Denote by $\mathcal{D}_{\mu}$ a data set consisting of $N$ samples for each of the $K$ classes – each sampled from the Gaussian distribution with centers $\mu_{1},\dots,\mu_{K}$, respectively, and with a common covariance matrix $\Sigma$. Let $\hat{\lambda}$ be the smallest $\lambda$-value with the minimal estimated expected generalization error, as determined by cross validation on $\mathcal{D}_{\mu}$. Denote by $\operatorname*{Err}_{\mu}(\lambda,\alpha)$ the generalization error of the model $\hat{\beta}(\lambda,\alpha)$ that has been estimated from the training set $\mathcal{D}_{\mu}$, by the sparse group lasso, for the given values of $\lambda$ and $\alpha$. Then let

where $\operatorname*{Err}_{\text{Bayes}}(\mu)$ is the Bayes rate. We are interested in trends in $Z_{\mu}$, as a function of $\alpha$, for different configurations of $\mu_{1},\dots,\mu_{K}$. To be specific, we will sample $\mu_{1},\dots,\mu_{K}$ from one of the following distributions:

• 1.

  A sparse model distribution, where the first $p_{a}$ entries of $\mu_{i}$ are i.i.d. with a distribution that is a mixture of the uniform distribution on $[-2,2]$ and the degenerate distribution at 0 with point probability $p_{0}$.

• 2.

  A dense model distribution, where the first $p_{a}$ entries of $\mu_{i}$ are i.i.d. Laplace distributed with location 0 and scale $b$.

The last $p_{b}$ entries are zero. We take $p_{a}=\lfloor 5/(1-p_{0})\rfloor$ throughout for the sparse model distribution. The within class covariance matrix $\Sigma$ is constructed using features from the cancer site data set. Let $\Sigma_{0}$ be the empirical covariance matrix of $p$ randomly chosen features. To avoid that the covariance matrix become singular we take

for $\delta\in(0,1)$.

The primary quantity of interest is

the expectation being over $\mu$ and the data set $\mathcal{D}_{\mu}$. We are also interested in how well we can estimate the non-zero patterns of the $\mu_{i}$’s. Consider this as $Kp$ two class classification problems, one for each parameter, where we predict the $\mu_{ij}$ to be non-zero if $\hat{\beta}_{ij}$ is non-zero, and $\mu_{ij}$ to be zero otherwise. We calculate the number of false positives, true positives, false negatives and true negatives. The positive predictive value (ppv) and the true positive rate (tpr) are of particular interest. The true positive rate measures how sensitive a given method is at discovering non-zero entries. The positive predictive value measures the precision with which the method is selecting the non-zero entries. We consider the following two quantities

In order to estimate the quantities (5) and (6) we sample $M$ configurations of $\mu$ from one of the above distributions. For each configuration we sample a training and a test data set of sizes $NK$ and $100K$, respectively. Using the training data set we fit the model $\hat{\beta}(\hat{\lambda},\alpha)$ and estimate $Z_{\mu}(\alpha)$ using the test data set. Estimates $\widehat{\operatorname*{err}}(\alpha)$, $\widehat{\operatorname*{tpr}}(\alpha)$ and $\widehat{\operatorname*{ppv}}(\alpha)$ are the corresponding averages over the $M$ configurations.

For this study we chose $M=100$, $N=15$, $K=25$, $p_{b}=50$, $\delta=0.25$ and the following three configuration distributions:

• 1.

  Thin configurations, where the centers are distributed according to the sparse model distribution with $p_{0}=0.95$, as defined above.

• 2.

  Sparse configurations, where the centers are distributed according to the sparse model distribution with $p_{0}=0.80$.

• 3.

  Dense configurations, where the centers are distributed according to the dense model distribution with scale $b=0.2$ and $p_{a}=25$.

In Figure 4 we see that for thin configurations the lasso has the lowest estimated error gab, along with the sparse group lasso with $\alpha=0.8$. For the sparse configurations the results indicate that the optimal choice of $\alpha$ is in the open interval $(0,1)$, but in this case all choices of $\alpha$ result in a comparable error gab. For the dense configurations the group lasso is among the methods with the lowest error gab.

In Figure 5 we plotted the true positive rate for the three configurations. Except for the thin configurations, the lasso is markedly less sensitive than the sparse group and group lasso methods. However, looking at Figure 6 we see that the sparse group and group lasso methods have a lower precision than the lasso, except for the dense configurations. We note that the group lasso has the worst precision, except for the dense configurations.

We provide two software solutions in relation to the current paper. An R package, msgl, with a relatively simple interface to our multinomial and logistic sparse group lasso regression routines. In addition, a C++ template library, sgl, is provided. The sgl template library gives access to the generic sparse group lasso routines. The R package relies on this library. The sgl template library relies on several external libraries. We use the Armadillo C++ library [14] as our primary linear algebra engine. Armadillo is a C++ template library using expression template techniques to optimize the performance of matrix expressions. Furthermore we utilize several Boost libraries [15]. Boost is a collection of free peer-reviewed C++ libraries, many of which are template libraries. For an introduction to these libraries see for example [16]. Use of multiple processors for cross validation and subsampling is supported through OpenMP [17].

The msgl R package is available from CRAN. The sgl library is available upon request.

We developed an algorithm for solving the sparse group lasso optimization problem with a general convex loss function. Furthermore, convergence of the algorithm was established in a general framework. This framework includes the sparse group lasso penalized negative-log-likelihood for the multinomial distribution, which is of primary interest for multiclass classification problems.

We implemented the algorithm as a C++ template library. An R package is available for the multinomial and the logistic regression loss functions. We presented applications to multiclass classification problems using three real data examples. The multinomial group lasso solution achieved optimal performance in all three examples in terms of estimated expected misclassification error. In one example some sparse group lasso solutions achieved comparable performance based on fewer features. If there is cost associated with the acquisition of each feature, the objective would be to minimize the cost while optimizing the classification performance. In general, the sparse group lasso solutions provide sparser solutions than the group lasso. These sparser solutions can be of interest for model selection purposes and for interpretation of the model.