Extreme Sparse Multinomial Logistic Regression: A Fast and Robust Framework for Hyperspectral Image Classification

Although the sparse multinomial logistic regression (SMLR) has provided a useful tool for sparse classification, it suffers from inefficacy in dealing with high dimensional features and manually set initial regressor values. This has significantly constrained its applications for hyperspectral image (HSI) classification. In order to tackle these two drawbacks, an extreme sparse multinomial logistic regression (ESMLR) is proposed for effective classification of HSI. First, the HSI dataset is projected to a new feature space with randomly generated weight and bias. Second, an optimization model is established by the Lagrange multiplier method and the dual principle to automatically determine a good initial regressor for SMLR via minimizing the training error and the regressor value. Furthermore, the extended multi-attribute profiles (EMAPs) are utilized for extracting both the spectral and spatial features. A combinational linear multiple features learning (MFL) method is proposed to further enhance the features extracted by ESMLR and EMAPs. Finally, the logistic regression via the variable splitting and the augmented Lagrangian (LORSAL) is adopted in the proposed framework for reducing the computational time. Experiments are conducted on two well-known HSI datasets, namely the Indian Pines dataset and the Pavia University dataset, which have shown the fast and robust performance of the proposed ESMLR framework.


Introduction
Although the rich spectral information available in a hyperspectral image (HSI) allows classifying among spectrally similar materials [1], supervised classification of HSI remains a challenging task, mainly due to the fact that unfavorable ratio between the limited number of training samples and the large number of spectral band [2].It may result in the Hughes phenomenon when the spectral bands increase and leads to poor classification results [3].To tackle such challenges, a number of state-ofthe-art techniques have been proposed, such as the support vector machine (SVM) [4], the multi-kernel classification [5], the extreme learning machine (ELM) [6] and the sparse multinomial logistic regression [2,[7][8][9][10] (SMLR).In addition, many approaches have also been proposed for dimensionality reduction and feature extraction [11][12], which include the principal component analysis and its variations [13][14][15][16] (PCA), the extended multi-attribute profiles [17][18] (EMAPs), the singular spectrum analysis [19][20][21][22] (SSA) and the segmented auto-encoder [12,22].Among these methods, the SMLR [23][24] has drawn a lot of attentions due to its good performance.
The SMLR has been proved to be robust and efficient under the Hughes phenomenon and is able to learn the class distributions in a Bayesian framework [25].Hence, it can provide a degree of plausibility for performing these classifications [26].Moreover, the logistic regression via the variable splitting and the augmented Lagrangian (LORSAL) has been proposed for dealing with large datasets and multiple classes efficiently.Since it can effectively learn a sparse regressor with a Laplacian prior distribution of the SMLR, the combination of SMLR and LORSAL is found to be one of the most effective methods for coping with the high dimensional data of HSI [26][27].
However, the existing SMLR framework suffers from some severe drawbacks.First, the SMLR with the original spectral data of the HSI as features is inefficient, thus it is necessary to find a better representation of the HSI data for improved classification accuracy.The second is the manually set initial value for the regressor, which may result in poor classification of HSI due to improper initial value used.Recently, some deep learning algorithms such as the convolutional neural network [28][29] (CNN) and the extreme learning machine (ELM) have drawn lots of attentions due to their good classification results for the HSI [30][31][32].However, CNN requires huge computational time and seem unrealistic.ELM is a generalized single layer feedforward neural networks (SLFNs), which characterizes fast implementation, strong generalization capability and a straightforward solution [6].The main goals of ELM are to minimize the output weights of the hidden layer and maintain the fast speed.According to the Bartlett's neural network generalization theory [33], the smaller norm of the weights will lead to better generalized performance.Hence, a feedforward neural network can reach a smaller training error [34].
For efficiency, the input weights and the bias between the input layer and the hidden layer of the ELM are randomly generated.It has been proved to be a fast and good data representation method [30][31][32].In fact, besides ELM, some other models, such as the liquid state machines [35] and the echo state networks [36][37] have also adopted the random weight selection technique with great success [38].Therefore, the problem of the poor data representation in the SMLR can be addressed using the random weight selection technique.Hence, in this paper, we propose the extreme sparse multinomial regression (ESMLR) for the classification of HSI.First, the data in the HSI will be represented by randomly generated weight and bias for SMLR, which also maintain the fast speed of the SMLR and improve the representation performance.Second, we set up an optimization model to minimize the training error of the regressor value, which is solved by using the Lagrange multiplier method and dual principle in order to automatically find a better initial regressor value for the SMLR (detailed in Section III).
In addition to spectral features, spatial information is also very important for the classification of the HSI.In the proposed ESMLR framework, the extended multi-attribute profile (EMAP) is used for feature extraction, as both morphological profiles (MPs) [39] and the attribute profiles (APs) [16][17] have been successfully employed for performing the spectral and spatial HSI classification.Moreover, the linear multiple feature learning (MFL) [7] is employed to maintain the fast speed and further improve the classification results.The MFL has been proposed for adaptively exploiting the information from both the derived linear and nonlinear features.As a result, it can potentially deal with the practical scenarios that different classes in the HSI datasets need different (either nonlinear or linear) strategies [7].According to the Li's works [7], the nonlinear feature such as the kernel feature contributes little to the HSI classification when the MFL is utilized.Moreover, it requires much more computational efforts for processing the nonlinear features.Therefore, a linear combination of the MFL which just utilizes the linear features of the HSI is proposed for the ESMLR.Hence, this operation can not only improve the classification results, but also maintain the fast speed of the ESMLR.
The main contributions of the proposed ESMLR framework in this paper can be highlighted as follows.First, the problem of the SMLR that uses the initial data of the HSI for performing the classification is addressed by randomly generating the input weights and bias of the input data, which will not only maintain the fast processing speed, but also improve the classification results of the HSI.Second, a new principle is introduced to automatically determine a suitable initial regressor value for SMLR to replace the manually settings.Third, the linear combination of the MFL that integrates the spectral and spatial information of HSI extracted by EMAPs followed by LORSAL is employed for the ESMLR for fast and robust data classification of HSI.
The remainder of this paper is structured as follows.Section II describes the background of the SMLR and the EMAPs.In Section III, the proposed ESMLR framework is presented.The experiment results and discussions are summarized in Section IV.Finally, Section V concludes this paper with some remarks and suggestions for the plausible futures.

EMAPs
The APs are obtained by applying attribute filters (AFs) to a gray-level image [16].The AFs are connected operators defined via a mathematical morphological mean for a gray level image to keep or to merge their connected components at different gray levels [39].Let  and  be an attribute thinning and an attribute thickening based on an arbitrary criterion   .Given an image   and a sequence of thresholds { 1 ,  2 , … ,   }, an AP can be obtained by applying a sequence of attribute thinning and attribute thickening operations as follows: AP(  ) = {   (  ),   −1 (  ), … ,   1 (  ),   ,   1 (  ), … ,   −1 (  ),   −1 (  )}. (1) Note that in (1) the AP is defined on each spectral band, hence the dimensionality of the APs will be very high when it is applied for the full spectral bands of the HSI [39].In [40], the principal component analysis (PCA) was suggested to solve this problem.Hence, the extended AP (EAP) is acquired by generating an AP on each of the first c PCs below [41].EAP = {( 1 ), ( 2 ), … , ( c )} Then, the EMAPs is defined as the composition of b different EAPs based on a set of b attributes { 1 ,  2 , . . .,   } as follows: Although a wide variety of attributes can be applied to the APs [42] for performing the HSI classification, in this paper we only consider the area attribute in order to maintain the fast speed whilst incorporating the spectral and spatial information.
Here, the code of the APs is from online http://www.lx.it.pt/~jun/.The threshold values of the area attribute were chosen as 100, 200, 500 and 1000.The first c PCs are determined to have the cumulative eigenvalues larger than 99% of the total value.

SMLR
Let t={1, …, M} be a set of M class labels.Denote S={1, …, n} as a set of integers indexing the n pixels of any HSI and x=( 1 , … ,   ) ∈  × be the HSI.Here, each pixel in the HSI is a d-dimensional feature vector and y = ( 1 , … ,   ) denote the labels of x.Let   = {( 1 ,  1 ), … , (  ,   )} be the training set.All the above parameters will be discussed in Section III.
First of all, the posterior class probabilities are modeled by the MLR [2], [7] as follows: where w = [ (1) , … ,  (−1) ]  ∈  −1×d denotes the regressors and h(  ) denotes the input feature.Here, the superscript 'T' denotes the transpose operator of a matrix.  is set to be 0 because the densities of ( 4) do not depend on the translation of the regressor   [7].The input features h can be linear or nonlinear.In the former case, we have: where  , is the j-th component of   .
If h(. ) is nonlinear, it can be formulated as follows: where φ(. ) is a nonlinear function.
According to [2], [7], the regressor w of the SMLR can be obtained by calculating the maximum a posteriori estimate as follows: Here, ℓ() is the logarithmic likelihood function given by: It is worth noting that log () is a prior over  which is irrelevant to the observation x.For controlling the complexity and the generalization capacity of the classifier,  is modeled as a random vector with the Laplacian density denoted as Here,  is the regularization parameter controlling the degree of sparsity [2].Hence, the solution of SMLR can be expressed as follows: The LORSAL algorithm is applied to SMLR to cope with the larger size problem of the HSI data.

From SMLR to ESMLR
The MLR can be modeled as follows [24], [43]: and where h(  ) is the input feature of the MLR and w = [ (1) , … ,  (−1) ]  ∈  (−1)×d denotes the regressors.  is set to be 0, as the densities of ( 10) and (11) do not depend on the translation of the regressor   [7].The input features h can be linear or nonlinear.
If the input feature is linear, then we have: where  , is the j-th component of   .
If h is nonlinear, it can be formulated as: where φ(. ) is a nonlinear function.
The initial regressor value can be used to find a better representation of the HSI for the ESMLR via determining the solution of the following optimization problem: where Y = [ 1 * , … ,   * ] ∈  (−1)× , and .
Here, if   belongs to the j-th class,  , * = 1.Otherwise,  , * = 0.In fact, the activation function ℎ() can be either linear or nonlinear, and L is the dimension of the feature space which we want to project;   ∈   and   ∈  1 are randomly generated.Actually, a wide range of feature mapping functions can be considered in our work which include but not limit to: (1) Linear function: ℎ(  ,   ,   ) =      +   ; (2) Sigmoid function: ℎ(  ,   ,   ) = From ( 14), it can be seen that the objective of the optimization is not only to reach a smaller training error, but also to reach a smaller value of the regressor w.According to the Bartlett's theory [33], this will help the proposed approach to achieve a good performance.From the optimization theory viewpoint [30]- [34], ( 14) can be reformulated as follows: where C is a regularization parameter and   is the training error for the samples   .
Based on the Karush Kuhn Tucker optimality conditions and the Lagrange multiplier method [44], we have: where  , is the Lagrange multiplier.
Then, the optimization condition can be expressed as follows: where Hence, the solution of the optimization defined in ( 15) can be analytically expressed as or From ( 20) and ( 21), it can be seen that the initial regressor value  0 is good for the ESMLR satisfying the optimization condition.Here, the random weight function h() can be used not only to find a better representation of the HSI data, but also to maintain the fast speed for the proposed framework.Based on the principles of the SMLR algorithm [23], the regressor  of the proposed ESMLR at the k-th iteration can be computed by the maximum a posterior estimate as follows: where ( −1 ) ∝ exp(−‖ −1 ‖ 1 ) and  is the regularization parameter for controlling the degree of sparsity [23].
The solution of ( 22) at the k-th iteration can be addressed by introducing the linear or nonlinear input features.That is, for the linear case: h(  ) = [ℎ( ,1 ), … , ℎ(  )]  ); for the nonlinear (kernel) case: h(  ) = [1,  1 (  ,  1 ), … ,  d (  ,   )]  , where  is a nonlinear function.Also, we have: Similar to K-SMLR [2], [26], the proposed ESMLR can also be extended to form a kernel-based ESMLR (K-ESMLR).The performance of the proposed ESMLR and K-ESMLR are evaluated in next section.In order to address the larger size problem of the HSI data, including large datasets and the number of classes, the LORSAL [45] algorithm is adopted.Moreover, the EMAPs are utilized for performing the efficient feature extraction and incorporating the spectral information and the spatial information.

ESMLR with a linear MFL
As mentioned above, the spectral information and the spatial information are integrated to further improve the performance of the proposed framework.It is well known that the kernel transform will increase the size of the input feature.As shown in [2], [7], the kernel transform may contribute slightly on the HSI classification accuracy when nonlinear features are utilized for the MFL.The kernel feature will also slow the speed of algorithms.Based on this perspective, a combinational linear MFL is proposed for improved HSI data classification whilst maintaining the low computational time of the proposed ESMLR.
Let ℎ  (  ) and ℎ  (  ) be the input spectral features of the raw/original HSI data and the spatial features extracted by the EMAPs, respectively.The input features of the proposed ESMLR can be expressed as follows: Then, (23) can be reformulated as: From (25), it can be seen that ( 23) and ( 25) have the same structure.Therefore, the LORSAL algorithm will be adopted in the proposed framework.Figure 1 shows the flowchart of the proposed spectral spatial ESMLR framework.

Experimental Section
In this section, the proposed ESMLR and K-ESMLR will be evaluated and relevant results are summarized in detail as follows.

The datasets
Two well-known HSI datasets are used in our experiments, which are detailed below.and water absorption bands.In total there are 42776 labelled samples in 9 classes within this dataset.

Compared Methods and parameter settings
The proposed ESMR framework are compared with the classical classifiers such as the K-SVM [34] (The codes of the K-SVM are obtained online from http://www.fst.umac.mo/en/staff/fstycz.html/), the SMLR and the K-SMLR [2], [7] (The codes of the SMLR and the K-SMLR are from online http://www.lx.it.pt/~jun/).All experiments are conducted in MATLAB R2015a and tested on a computer i7 with 3.40GHz CPU and 8.0G RAM.All data are normalized via the unit max method, i.e. each data of a HSI is divided by the largest value of the whole dataset.
For all kernel-based/nonlinear methods, the Gaussian radial basis function (RBF) kernel is used.For the parameter σ of the RBF in the K-SMLR and the K-ESMLR, it is set to be 0.85 for the Indian Pines dataset and 0.35 for the Pavia University dataset as suggested by Sun et al [26].The LIBSVM toolbox of the MATLAB R2015a is used for the implementation of the K-SVM approach [46].The parameter of the K-SVM were chosen according to [34].For the cost parameter in (20) or ( 21), C = 2 a is chosen where a is in the range {1, 2, …, 20}.The regularization parameter in ( 22) is set to  = 2  .The total number of dimension of the new feature space L is chosen in the range {50, 100, …, 1450, 1500}.If there is no special emphasis required, the dimension of the new feature space in the proposed ESMLR is set to be L=300 for spectral information only and the combined spectral-spatial information (EMAPs).Similarly, L=500 is chosen for the situation that utilize spectral and spatial information (linear MFL).For other parameters in the SMLR, K-SMLR, ESMLR and K-ESMLR, details will be discussed in the subsections below.All experiments are repeated 10 times with the average classification results reported for comparison.

Discussions on the robustness of the ESMLR framework
In the following experiments, the robustness of the proposed framework is evaluated.For the Indian Pines dataset and the Pavia University dataset, in total 515 and 3921 samples are randomly selected for training, respectively.The remaining samples are used for testing based on the overall accuracy (OA) of classification.For the two datasets, the number of samples used for training and testing from each class are summarized in Table 1 and Table 2, respectively, along with the classification results under different experiments settings as detailed below.
Experiment 1: In this experiment, the proposed ESMLR approach is evaluated in three different situations, i.e. using only the spectral information as features, using both spectral and spatial information yet in combination with EMAPs and the linear MFL, respectively.For the Indian Pines and the Pavia University datasets, the results are shown in Figure 2 and Figure 3, respectively.As seen, the proposed framework shows the good performances in all the situations when L>150.The fusion of spectral and spatial information can successfully improve the OA, where the combination of ESMLR and MFL slightly outperforms ESMLR and EMAPs.
Experiment 2: In this experiment, the impact of the parameter C ( C = 2 a ) in the proposed K-ESMLR under the aforementioned three different situations are evaluated.As shown in Figure 4, the proposed ESMLR achieves a very good performance when C is larger than 0, where the classification results are very stable even though a or C is significantly changed.
This again demonstrates the robustness of the proposed framework.

Experiment 3:
The impact of the sparse parameter b ( = 2  ) on ( 22) are evaluated in this experiment.As shown in Figure 5, the proposed K-ESMLR achieves better classification results compared with the SMLR when only spectral information is utilized, especially for the Pavia University dataset.Also it seems K-ESMLR slightly outperforms K-SMLR.This has demonstrated that the proposed framework achieves a better performance compared with the conventional SMLR framework for both linear and nonlinear (kernel) cases.When EMAPs are applied to extract both the spectral and spatial information, the proposed framework also achieves better classification results compared with SMLR.When combining our ESMLR framework with the proposed combinational linear MFL, it also outperforms SMLR.Note that K-SMLR and K-ESMLR cannot be combined with linear MFL, hence the results are not shown.In summary, thanks to the proposed improvements, the proposed ESMLR/K-ESMLR framework has outperformed conventional SMLR and K-SMLR for effective classification of HSI. (

Discussions on classification results and the running time of different algorithms
In this subsection, the classification results and the running (executed) time based on the proposed classifiers are compared with other state-of-the-art approaches.For Indian Pines and the Pavia University datasets, the results are summarized in Table 1 and Table 2, respectively.It is worth noting that all the classification results were based on the corresponding best parameters.3 and Table 4, respectively.As seen in Table 3 and Table 4, the proposed ESMLR framework improves the classification accuracy of SMLR dramatically even for a small number of training samples.When both spectral and spatial information are utilized, the proposed framework outperforms K-SVM.As K-SVM requires much more computational time in comparison to the proposed framework, we can conclude that the proposed ESMLR framework provides a fast and robust solution for the classification of HSI.

Classification results with different numbers of training samples
In addition, Figure 6 and Figure 7 show the classification results from different classifiers for the Indian Pines dataset and the Pavia University dataset, respectively.For each class, the number of training samples is set to 40.Again, this has clearly shown the superior performance of the proposed approach in effective classification of HSI.

Conclusions
In this paper, we propose a new ESMLR framework to solve the two main drawbacks SMLR for the effective classification of the HSI.By combining linear MFL for incorporating the spectral and spatial information of HSI, the classification accuracy has been successfully improved.Compared with conventional SMLR method, the proposed ESMLR framework has yielded better classification results with a comparable computational time.In comparison to K-SVM, ESMLR requires much less computation time and can better exploit the combination of spatial and spectral information with different labeled numbers of training samples.Furthermore, the proposed approach consistently achieves higher classification accuracy even under a small number of training samples.
The future works will focus on the optimization of the required computational time for K-ESMLR by using sparse representation, and further improvement of the classification accuracy by resorting the ideal regularized composite kernel [47].

Figure 1 .
Figure 1.The flowchart of the proposed ESMLR framework.

1 )
. The Indian Pines dataset: The HSI image was acquired by the AVRIS sensor in 1992.The image contains 145×145 pixels and 200 spectral bands after removing 20 bands influenced by the atmospheric affection.There are 10366 labelled samples in 16 classes within the HSI dataset.2).The Pavia University dataset: The system was built by the University of Pavia of Italy in 2001.The Pavia University dataset was acquired by the ROSIS instrument.The image contains 610×340 pixels and 103 bands after discarding 12 noisy

Figure 2 .Figure 3 .
Figure 2. The robustness performance of the proposed framework based on the Indian Pines dataset: (a) The proposed ESMLR with spectral information (200 features and b=-7); (b) The proposed ESMLR with spectral and spatial information (EMAPSs) (36 features and b=-11); (c) The proposed ESMLR with spectral and spatial information (proposed linear MFL) (236 features and b=-10).

Figure 4 .Figure 5 .
Figure 4.The robustness of the proposed framework under different values of cost C: (a) The proposed ESMLR with the Indian Pines dataset (using the original HSI dataset, i.e. 200 features, b=-15) and the K-ESMLR with the EMAPs (36 features, b=-17); (b) The proposed K-ESMLR with the Pavia University dataset (using the original HSI dataset, i.e. 103 features, b=-10) and the K-ESMLR with the EMAPs (36 features, b=-12).
In this section, we evaluate the robustness of the proposed framework with different numbers of training samples.We vary the number of training samples Q randomly selected from each class, where we have Q=5, 10, 15, 20, 25, 30, 35 and 40 in our experiments.If Q becomes more than 50% of the total samples within a class, only 50% of samples within that class is used for training.For the Indian Pines and the Pavia University datasets, relevant results are summarized in Table

Table 3 .
Classification accuracy (%) with different numbers of labeled samples in Indian Pines dataset (Best result of each row is marked in bold type).

Table 4 .
Classification accuracy with different numbers of labeled samples in Pavia University dataset (Best result of each row is marked in bold type).