Feature Extraction for Hyperspectral Imagery: The Evolution from Shallow to Deep
Abstract
The final version of the paper can be found in IEEE Geoscience and Remote Sensing Magazine.
Hyperspectral images provide detailed spectral information through hundreds of (narrow) spectral channels (also known as dimensionality or bands) with continuous spectral information that can accurately classify diverse materials of interest. The increased dimensionality of such data makes it possible to significantly improve data information content but provides a challenge to the conventional techniques (the socalled curse of dimensionality) for accurate analysis of hyperspectral images. Feature extraction, as a vibrant field of research in the hyperspectral community, evolved through decades of research to address this issue and extract informative features suitable for data representation and classification. The advances in feature extraction have been inspired by two fields of research, including the popularization of image and signal processing as well as machine (deep) learning, leading to two types of feature extraction approaches named shallow and deep techniques. This article outlines the advances in feature extraction approaches for hyperspectral imagery by providing a technical overview of the stateoftheart techniques, providing useful entry points for researchers at different levels, including students, researchers, and senior researchers, willing to explore novel investigations on this challenging topic. In more detail, this paper provides a bird’s eye view over shallow (both supervised and unsupervised) and deep feature extraction approaches specifically dedicated to the topic of hyperspectral feature extraction and its application on hyperspectral image classification. Additionally, this paper compares 15 advanced techniques with an emphasis on their methodological foundations in terms of classification accuracies. Furthermore, to push this vibrant field of research forward an impressive amount of codes and libraries is shared at: https://github.com/BehnoodRasti/HyFTechHyperspectralShallowDeepFeatureExtractionToolbox.
I Introduction
Hyperspectral imaging technology provides detailed spectral information by sampling the reflective portion of the electromagnetic spectrum covering a wide range from the visible region (0.40.7 Î¼m) to the shortwave infrared (SWIR) region (almost 2.4 Î¼m). Hyperspecral sensors can also characterize the emissive properties of objects by acquiring data in the range of the midwave and longwave infrared regions, in hundreds of narrow contiguous spectral channels.
Detailed spectral information provided by hyperspectral sensors present both challenges and opportunities. For instance, hyperspectral images can be used to differentiate between different classes of interest with slightly different spectral characteristics [26]. However, most of the commonly used methods utilized for the analysis of gray scale, color, or multispectral images cannot be extended to analyse hyperspectral images due to several reasons, which will be detailed later in Section I.A.
The limited availability of training samples (which is a common issue in remote sensing) dramatically impacts the performances of supervised classification approaches due to the high dimensionality of hyperspectral images, which poses a problem for designing robust statistical estimations. Feature extraction can be used as a solution to address the aforementioned issue, which can be described as finding a set of vectors that represents an observation while reducing the dimensionality by transforming the input data linearly or nonlinearly to another domain and extract informative features in the new domain. The use of feature extraction techniques can be advantageous for a number of reasons which is illustrated in Fig. 1 and described in the following subsections.
Ia Unique Geometrical and Statistical Properties of High Dimensional Data and the Need for Feature Extraction
Several studies (e.g., [98, 112, 51]) have demonstrated the unique geometrical, statistical, and asymptotical properties of highdimensional data compared to RGB and multispectral images. These properties, which have been shown through some experimental and theoretical examples, clearly justify the reasons why most of the analytical approaches that are developed for RGB and multispectral images are not applicable to hyperspectral images [6]. Among those experimental examples, we can recall that (1) as dimensionality increases, the volume of a hypercube concentrates in corners, or (2) as dimensionality increases, the volume of a hypersphere concentrates in an outside shell. With respect to these examples, the following conclusions have been drawn:

A highdimensional feature space is almost empty, which indicates that multivariate data in ( represents the number of bands, spectral channels, or dimensions) can usually be represented in a lower dimensional space (referred to as subspace) without loosing considerable information in terms of class separability [6].

Since the highdimensional feature space is almost empty (i.e., Gaussian distributed data have a tendency to concentrate in the tails while uniformly distributed data have a tendency to be concentrated in the corners), the density estimation of hyperspectral data for both Gaussian and uniform distributions become extremely challenging.
Fukunaga [21] claimed that there is a relation between the type of the classifier, the required number of training samples, and the number of input dimensions. As reported in [21], the required number of training samples is linearly related to the dimensionality for linear classifiers and to the square of the dimensionality for quadratic classifiers (e.g., the Gaussian maximum likelihood classifier [21]), and for nonparametric classifiers, the number of required samples exponentially increases as the dimensionality increases. Landgrebe showed a groundbreaking fact that too many spectral bands might have negative impacts in terms of expected classification performance [63]. When dimensionality increases, with a constant and limited number of training samples, a higher amount of statistics must be estimated. Thus, the accuracy of the statistical estimation decreases although higher spectral dimensions increase the separability between the classes. This leads to a decrease in classification accuracies beyond an unknown number of bands. These problems are related to the curse of dimensionality, also known as Hughes phenomenon [46]. This finding was against the general understanding of hyperspectral data where it was wrongly believed that full dimensionality is always better than subspace in terms of classification accuracies.
The unique characteristics of highdimensional data, as discussed above, have a pronounced impact on the performances of supervised classifiers [87], as they demand an adequate number of training samples, which is almost impossible to obtain in practice since the collection of such training samples is either expensive or time demanding. To address this issue, feature extractionbased dimensionality reduction is found to be effective.
IB Storage Systems and Processing Times
We are now in the era of massive data acquisitions. Statistics demonstrate that the cumulative volume of existing big data has been tremendously increased from 4.4 ZB to 44 ZB from 2013 to 2020
Feature extractionbased dimensionality reduction helps in data compression, which leads to the reduced storage space, faster transmission time, removing redundant features, reducing the storage space required, and fasten the required time for performing the same computations.
IC An EverGrowing Relation between Machine learning and Feature Extraction
Fig. 2 illustrates the basic idea of a machine learning approach, which consists of feature extraction and classification. In machine learning, users are requested to provide some guidelines for the machine (algorithm). These guidelines are usually provided by applying handcrafted feature extraction approaches to provide informative features for the subsequent classifier. At the very beginning, each image pixel is regarded as a pattern, and its spectrum (i.e., a vector of different values of a pixel in different spectral channels) is considered as the initial set of features. This set of features, which are also known as spectral features, suffers from two important downsides: (1) These features are often redundant and (2) they do not consider the spatial dependencies of the adjacent pixels.
To address the first issue, a feature reduction (through feature extraction or selection) step can be applied, aiming at reducing the dimensionality of the input data (from dimensions in the original data to dimensions in a new feature space ). This step can also be known as spectral feature extraction, which tries to preserve the key spectral information of the data by reducing the dimensionality and maximizing separability between classes. It is interesting to note that the second issue can also be addressed using feature extraction approaches. Please note that here the feature extraction step, which is also known as spatial feature extraction, is not aiming at reducing the dimensionality and instead is trying to model (extract) spatial contextual information suitable for the subsequent classification or object detection step and usually leads to an increase in the number of features. The simultaneous use of spectral and spatial features for hyperspectral data classification were studied in numerous works such as [6, 25, 124, 20].
Deep learning (as shown in Fig. 2), which is regarded as a subset of machine learning, tries to automatize the building blocks of the machine learning approaches (i.e., feature extraction and classification) by developing an endtoend framework, which takes the input, performs automatic feature extraction and classification by considering the unique nature of the input data (instead of those handcrafted feature extraction designs in machine learning), and outputs classification maps. It turns out that if an adequate amount of training data is supplied, deep learning approaches can outperform any other shallow machine learning approaches in terms of accuracies. Here, a question arises: Due to the fact that the number of available training data is often limited in the remote sensing community, would advanced deep learningbased approaches outperform their shallow alternatives in terms of accuracies? This question will be addressed in this paper.
Based on the abovementioned descriptions, feature extraction is a key step in both machine learning and deep learning, whose concept has been evolved significantly through time from unsupervised to (semi)supervised, from spectral or spatial to spectralandspatial, from manual to automatic, from handcrafted to endtoend, and from shallow to deep.
ID Contributions
This article provides a detailed and organized overview of hyperspectral feature extraction techniques, categorized into two general sections: shallow feature extraction techniques (further categorized into supervised and unsupervised) and deep feature extraction techniques. Each section provides a critical overview of the state of the art that is mainly rooted in the signal and image processing, statistical inference, and machine (deep) learning fields. Then, a few representative and advanced feature extraction approaches are chosen from each of the abovementioned categories for further analysis and comparisons (mostly in terms of usefulness for classification). This article will, therefore, contribute to answering the following questions: When it comes to hyperspectral data in Earth observation, are deep learningbased feature extraction approaches better alternatives than their conventional (yet advanced) shallow feature extraction techniques? Which factors should be considered to design robust shallow and deep feature extraction techniques? In addition, to further promote this field of research, this paper is accompanied with a significant amount of codes and libraries for hyperspectral feature extraction, which is made publicly available at https://github.com/BehnoodRasti/HyFTechHyperspectralShallowDeepFeatureExtractionToolbox. Finally, several possible future directions are highlighted. To make the contribution of this paper clearer compared to the existing papers in the literature, here we provide a brief discussion. The work of [68] is dedicated to the evolution of discriminant analysisbased feature extraction models, which is a specific type of dimensionality reduction approaches. The work of [49] reviewed feature extraction and data mining works, which had been published mostly until 2012. Since 2012, however, many deep and shallow feature extraction approaches have been developed, which are critically reviewed and compared against each other in this work. The work of [104] focuses only on feature selection approaches while our proposed paper covers feature extraction techniques, and therefore, they complement each other.
Ii Datasets and Notations
Iia Datasets
Indian Pines 2010
This dataset (Fig. 3) is a very high resolution hyperspectral image (VHR HS) acquired by the ProSpecTIR system over near Purdue University, Indiana, between 2425th of May 2010. In this paper, we use a subset of 445750 pixels with 360 spectral bands. The dataset has the spatial resolution of 2 m and spectral width of 5 nm. The dataset contains 16 land cover classes shown in Fig. 3. The training and test sets used in this study have also been demonstrated in Fig. 3. Table I gives the number of samples, including training and test samples used in the experimental section.
Class No.  Class Name  Training Samples  Test Samples  Samples 

1  Cornhigh  726  2661  3387 
2  Cornmid  465  1275  1740 
3  Cornlow  66  290  356 
4  Soybeanhigh  324  1041  1365 
5  Soybeanmid  2548  35317  37865 
6  Soybeanlow  1428  27782  29210 
7  Residues  368  5427  5795 
8  Wheat  182  3205  3387 
9  Hay  1938  48107  50045 
10  Grass/Pasture  496  5048  5544 
11  Cover crop 1  400  2346  2746 
12  Cover crop 2  176  1988  2164 
13  Woodlands  1640  46919  48559 
14  Highway  105  4758  4863 
15  Local road  52  450  502 
16  Buildings  40  506  546 
Total  10954  187120  198074 
Houston University 2012
This dataset was acquired on June 23, 2012, by the Compact Airborne Spectrographic Imager (CASI) over the campus of the University of Houston and the neighboring urban area. The average height of the sensor was 5500ft. The data contain 349 1905 pixels with the spatial resolution of 2.5 m and 144 spectral bands ranging 0.381.05 . The dataset includes 15 classes of interests shown in Fig. 4. A color composite representation of the data and the corresponding training and test samples used in this study are shown in Fig. 4. The number of training and test samples for different classes of interests used in the experiments are given in Table II.
Class No.  Class Name  Training Samples  Test Samples  Samples 

1  Grass Healthy  198  1053  1251 
2  Grass Stressed  190  1064  1254 
3  Grass Synthetic  192  505  697 
4  Tree  188  1056  1244 
5  Soil  186  1056  1242 
6  Water  182  143  325 
7  Residential  196  1072  1268 
8  Commercial  191  1053  1244 
9  Road  193  1059  1252 
10  Highway  191  1036  1227 
11  Railway  181  1054  1235 
12  Parking Lot 1  192  1041  1233 
13  Parking Lot 2  184  285  469 
14  Tennis Court  181  247  428 
15  Running Track  187  473  660 
Total  2832  12197  15029 
Houston University 2017
This dataset was acquired on Feb. 16, 2017 by the hyperspectral imager CASI 1500 over the area of the University of Houston. The data covers the spectral range 3801050 nm with 48 bands with the ground sampling distance of 1 m. It contains 601 2384 pixels and 20 land cover classes of interest shown in Fig.5. The VHR RGB image is downsampled and shown in Fig.5 together with the corresponding training and test samples used in this study. The number of training and test samples for different classes of interest used in the experiments are given in Table III.
Class No.  Class Name  Training  Test  Sample 

1  Healthy grass  1458  8341  9799 
2  Stressed grass  4316  28186  32502 
3  Synthetic grass  331  353  684 
4  Evergreen Trees  2005  11583  13588 
5  Deciduous Trees  676  4372  5048 
6  Soil  1757  2759  4516 
7  Water  147  119  266 
8  Residential  3809  35953  39762 
9  Commercial  2789  220895  223684 
10  Road  3188  42622  45810 
11  Sidewalk  2699  31303  34002 
12  Crosswalk  225  1291  1516 
13  Major Thoroughfares  5193  41165  46358 
14  Highway  700  9149  9849 
15  Railway  1224  5713  6937 
16  Paved Parking Lot  1179  10296  11475 
17  Gravel Parking Lot  127  22  149 
18  Cars  848  5730  6578 
19  Trains  493  4872  5365 
20  Seats  1313  5511  6824 
Total  34477  470235  504712 
IiB Notations
The observed HSI is denoted by where and are the number of spectral bands and pixels in each band, respectively. indicates the dimension of the feature space (the subspace). where denotes the matrix which contains the training samples. denotes the vector which contains the class labels where and denotes the number of classes. is the identity matrix and is the estimate of matrix . The Frobenius norm is denoted by . denotes the trace of matrix . The definitions of the symbols used in the paper are given in Table IV.
Symbols  Definition 

the th entry of the vector .  
the th entry of the matrix .  
the th column of the matrix .  
the th row of the matrix .  
norm of the vector i.e. the number of nonzero entries.  
norm of the vector , obtained by .  
norm of the vector , obtained by .  
norm of the matrix , obtained by .  
Frobeniusnorm of the matrix , obtained by .  
the estimate of the matrix .  
the trace of the matrix .  
total variation norm of the matrix is obtained by 
Iii Shallow Feature Extraction Techniques
Iiia Unsupervised Feature Extraction Techniques
Unsupervised feature extraction (UFE) often refers to the FE techniques which do not incorporate the knowledge of the ground (ground reference or labeled samples) to extract features. UFE techniques often rely on intrinsic characteristic of the HSI data such as geometric, spatial or spectral information to extract the features. Arguably, the main advantage of UFE compared with the other FE techniques is the lack of need for the training samples, which is of great importance in the case of remote sensing datasets. In this paper, four major UFE groups widelyused for HSI analysis are studied, which are categorized in the following subsections. Fig. 6 illustrates the graphical abstracts of those groups. Before explaining the UFE techniques in more details, we briefly refer to three groups of widelyused FE techniques which could also be assumed as UFE, however, will not be studied in details in this paper due to their specific applications. The first group includes a range of approaches such as normalized differential vegetation index (NDVI) and normalized differential water index (NDWI) which often rely on the knowledge of the characteristics of the sensors. The second group includes unmixing techniques, which could also be assumed as UFE techniques. They often exploit optimization techniques to show the fractions of materials existing in pixels based on some assumptions on the spectral signatures of the materials. Therefore, the final features extracted represent different materials in the scene at the subpixel level [7]. The third group includes an impressive number of approaches based on mathematical morphology, which hierarchically extract spatial and contextual information from the input image and usually leads to a significant increase in the number of features [24].
Conventional Data Projection/Transformation Techniques
Numerous UFE techniques fall into this category. The conventional techniques categorized in this group are often designed to linearly project or transform the data, , in a lower dimensional feature space (also called subspace) exploiting different nonlocal intrinsic characteristics of the hyperspectral dataset. The transformation can be given by
(1) 
where is the projected data in the lower dimensional space and is the transformation matrix or the bases for the subspace. Arguably, principal component analysis (PCA) [54] can be considered as the most conventional UFE technique which has been widely used for hyperspectral analysis [99]. PCA captures the maximum variance of the signal by projecting the signal on the eigenvectors of the covariance matrix () using
(2) 
A widely used HSI UFE technique is the maximum noise fraction (MNF) [29] or noise adjusted principal components (NAPC) [65] which seek a projection in which the signal to noise ratio (SNR) is maximized. MNF uses the following optimization
(3) 
where is the noise covariance matrix. Another conventional technique is independent component analysis (ICA) [47]. ICA assumes a linear mixture model of the nonGaussian independent source signals and the mixing matrix which both are simultaneously estimated and therefore ICA is referred to as blind source separation. ICA has been also widely used for hyperspectral image analysis [110].
To cope with the nonlinearity of the HSI data, the kernel (nonlinear) versions of the aforementioned techniques i.e., kernel MNF [84], kernel ICA (KICA) [1], and kernel PCA (KPCA) [97] have been also proposed. By using the kernel trick the data are projected into a feature space where the inner product are defined using a kernel function. KICA and KPCA have been used as UFE techniques for change detection and classification in [80] and [19], respectively. In [10], discrete wavelet transformation (DWT) has been used for hyperspectral feature extraction. Since, DWT does not reduce the dimension, in [10], linear discriminant analysis (LDA) has been exploited to reduce the dimension.
Band Clustering/Splitting and Mergingbased Techniques
The top right subfigure from Fig. 6 shows the basic steps of band clustering and mergingbased feature extraction methods. As shown in this figure, the core idea behind this group of methods is to split the spectral bands into several groups in which the spectral bands have very high correlation. Hence, the proposed techniques often use similarity and dissimilarity criteria to split the spectral bands into several nonoverlapping groups. By selecting or fusing the bands of each group, some representative bands or features of different groups are obtained. Furthermore, followed by the merging step, some band filtering and processing operations can be also used to further improve the discrimination of the resulting features. This group of techniques is often computationally cheap, and thus has been widely used in real applications. On the other hand, spectral information is often neglected by the methods in this category.
For the band clustering and merging two algorithms are proposed in [62]. The first algorithm selects discriminative bases by considering all the classes simultaneously, however, the second algorithm selects the best bases for a pair of classes at a time. In [81], a hierarchical clustering algorithm was introduced to split and cluster the hyperspectral bands where the representative band for each cluster is selected based on both a mutual information criterion and a divergencebased criterion. Another band clustering technique was proposed in [12] where the splitting is done by minimizing a mutual information criterion applied on averaged bands iteratively. Iterative algorithms were proposed in [88], for both splitting and merging the bands. The splitting procedure is done using the Pearson correlation coefficient between adjacent bands. Then, the merging is applied by averaging over the splited bands.
Besides splitting/clustering and merging hyperspectral bands, another operation is to further improve the feature discrimination by band filtering or processing. For example, a hyperspectral feature extraction using image fusion and recursive filtering was given in [55] where the adjacent bands are fused by averaging and then recursive filtering was used for extracting spatial information. In [57], the intrinsic image decomposition is applied for processing the merged bands, which can effectively remove information that is not related to the material of different objects. After that, multiple improved versions of intrinsic decompositionbased band processing methods are developed [52, 53]. In [17], a relative total variationbased structure extraction method is applied for band processing, so as to construct multiscale structural features which are robust to image noise.
LowRank Reconstructionbased Techniques
Lowrank reconstructionbased feature extraction techniques proposed by Rasti et al. [94, 89, 90, 92] are based on finding an orthogonal subspace by minimizing a constrained cost function. They exploit lowrank models and reconstructionbased optimization frameworks to extract features. The optimization frameworks take into account the prior knowledge of the data using different types of penalties. Due to the noise assumption in the lowrank model used, this group of FE techniques is robust to noise. They are often computationally expensive compared to groups 1 and 2 due to the iterative algorithms used to solve the (nonconvex) optimization problem.
Waveletbased sparse reduced rank regression (WSRRR) [89] applies the sparsity prior on the wavelet coefficients considering that the projected data on wavelet bases are sparse. WSRRR uses the model
(4) 
where represents 2D wavelet bases, is the observed HSI, contains the orthogonal subspace bases, and is the noise and model error. WSRRR simultaneously estimates the lowrank projection matrix and the wavelet coefficients which minimizes
(5) 
Note that the extracted features are given by .
To capture the spatial (neighboring) information, orthogonal total variation component analysis (OTVCA) has been proposed in [90] where the HSI is modeled as
(6) 
where matrix contains the unknown features. OTVCA assumes that the hyperspectral features are spatially piecewise smooth and, therefore, exploits the total variation (TV) penalty and simultaneously estimates and using
(7) 
where
and and are the matrix operators to calculate the first order vertical and horizontal differences, respectively, of a vectorized image. Recently, sparse and smooth lowrank analysis (SSLRA) was proposed in [92] which models the HSI based on a combination of sparse and smooth features
(8) 
where and contain smooth and sparse features, respectively. SSLRA extracts simultaneously the sparse features, , and the smooth ones, , by taking into account both sparsity and TV penalties
(9) 
Graph Embedding and/or Manifold Learning Techniques
Considering the nonlinear characteristic of HSIs, this group of FE techniques aims to capture the data manifold through local geometric structure of neighboring pixels in the feature space. Fig. 6 (bottom right) demonstrates the concept of manifold learning FE techniques applied on the Swiss roll dataset. The pink line in the left figure shows the Euclidean distance between two data points in 3D space. It is clear that this line is not an effective metric to measure the similarity of the two points selected in the Swiss roll dataset. On the other hand, after unfolding the dataset which can be represented in 2D space, in the right figure, the Euclidean distance between two data points shown by the pink line is a better representative of the similarity of the two point in the dataset. The FE techniques categorized in this group are indeed designed to capture such a manifold while representing the data in a lower dimensional feature space.
Graph embedding or manifold learning FE techniques often include three main steps, 1) neighborhood pixel selection 2) weight selection, and 3) embedding construction. Isometric mapping (ISOMAP) [109, 2] is a global geometric nonlinear feature extraction. ISOMAP searches for geodesic distances between data points. It includes three main steps; 1) Constructing a neighborhood graph of the data points. 2) Computing the shortest path distances between all data points in the neighborhood graph. 3) constructing the lower dimensional embedding vectors which preserves the path distances in the neighborhood graph. Locally linear embedding (LLE) [95], Laplacian eigenmaps (LE) [4], and locality preserving projection (LPP) [35] are also geometric nonlinear feature extraction based on graph embedding. LLE constructs the embedding graph in three steps; 1) Selecting the neighbors for data points using the nearest neighbors. 2) Compute the weights that linearly reconstruct the data points using their neighbors by minimizing the following constrained leastsquares
(10) 
where contains the neighborhood pixels selected for . We should note that the constrained weights estimated from (10) for every data point are invariant to rotations, rescalings, and translations of that data point and its neighbors and therefore they characterize the intrinsic geometric properties of each neighborhood. 3) Constructing the lower dimensional embedding vectors by minimizing
(11) 
We should note that the reconstruction weights are fixed in minimization (IIIA4) and therefore the intrinsic geometric properties of the data with dimension are invariant to such a transformation into a lower dimension .
In [118], a general framework for graph embedding is given by
(12) 
or equivalently
(13) 
where denotes the Laplacian matrix of the undirected weighted graph (where is the vertex set and is the similarity matrix), is a diagonal matrix where its entries are given by
(14) 
The diagonal matrix is for the scale normalization and might also be the Laplacian matrix of a penalty graph such as . We should note that the vertices of and (i.e., ) are the same while the similarity matrix () corresponds to the similarity characteristics suppressed in the lower dimensional feature space (, see [118]). LLE can be reformulated using graph embedding mentioned above with similarity matrix if otherwise and [118]. ISOMAP, LE, and LPP can be also formulated using graph embedding [118]. In the viewpoint of graph embedding the main differences between these FE techniques are the selection of the matrices and . For instance, LE and LPP use the Gaussian function with the standard deviation to choose the similarity matrix as
(15) 
We should note that the techniques categorized in this group are assumed as supervised FE methods when they are applied only on the training samples. This is common in the case of HSI due to the large volume of the image which makes the algorithm computationally very expensive. In the following section we will discuss how the ground reference (training samples) can be used to construct the edge matrix, W, and therefore those techniques are considered as SFE.
IiiB Supervised Feature Extraction Techniques
Unlike unsupervised FE techniques that rely on modeling various prior assumptions of hyperspectral data, supervised methods are capable of extracting classseparable features more effectively, owing to the use of label information. Over the past few decades, some seminal models have been widely developed and applied to perform supervised feature extraction on HSIs, which can be roughly categorized into two streams: Subspace learning (SL)based and band selection (BS)based approaches.
Different from the handcrafted features [114], the SLbased approaches learn to extract the lowdimensional representation from the data by formulating different supervised rules in view of label information. There are some typical methods in SL, including linear discriminant analysis (LDA) [18], matrix discriminant analysis (MDA) [33], decision boundary FE (DBFE) [64], etc. While the latter, BS, which aims at screening out the representative and informative spectral bands, is unfolded with mutual informationbased BS [31], rough set, and fuzzy Cmeans [101], to name a few. To further enhance the class separability, a large number of extended methods have been successfully proposed in recent years, which are subspace LDA (SLDA) [119], regularized LDA [3], local fisher’s discriminant analysis (LFDA) [103], feature space discriminative analysis (FSDA) [48], roughsetbased BS [86], and FE with local spatial modeling [11].
According to the powerful learning ability of SL methods compared to that of BSbased strategies, we rather focus on reviewing the SLrelated FE techniques, in which two main streams – discriminant analysis FE (DAFE) and regressioninduced representation learning (RIRL) – are emphatically investigated and compared by clarifying their similarities and differences as well as pros and cons, as briefly illustrated in Fig. 7.
Discriminant Analysis Feature Extraction (DAFE)
Generally speaking, DAFE seeks to find an optimal projection or transformation matrix ( is the dimension of the tobeestimated subspace) by optimizing certain classrelevant separation criterion associated with the label information. In this process, the estimated subspace that consists of a series of vector can be obtained by projecting the samples onto a decision boundary, which can be generally expressed as . Each vector in can be collected by . Depending on the different types of label embedding, DAFE can be subdivided into LDA and its variants, graphbased discriminant analysis (GDA) and its extensions, and kernelized discriminant analysis (KDA).
LDA and Its Variants: The traditional LDA linearly transforms the original data into a discriminative subspace by maximizing the Fisher’s ratio in the form of generalized Rayleigh quotient, that is, minimizing the intraclass scatter and maximizing interclass scatter simultaneously. Given a pairwise training set , the objective function of multiclass LDA to estimate the linear projection matrix can be written as follows:
(16) 
where and are defined as the withinclass scatter matrix and the betweenclass scatter matrix, respectively. With the constraint of , the optimization problem in (16) can be equivalently converted to one of by introducing the Lagrange multiplier . The closeform solution to the simplified optimization problem can be deduced by a generalized eigenvalues decomposition (GED).
Due to the sensitivity to complex highdimensional noises caused by the environmental and instrumental factors and the availability of labeled samples, the original LDA inevitably suffers from an illposed statistical degradation, especially in the case of smallscale samples. The degraded reasons mainly lie in the singularity of the two scatter metrics ( and ), thereby easily leading to the overfitting problem. To improve the stability and generalization, the regularized LDA was proposed by additionally adding a norm constraint on parameterized by as . By replacing in (16) with the regularized , the solution in the regularized LDA can be still obtained by the GED solver.
Considering the local neighborhood relations between samples in the process of model learning, LFDA breaks through the bottleneck of those LDAbased methods by assuming that the data are distributed in the nonlinear manifolds rather than a homogeneous Gaussian space. For this purpose, LFDA is capable of effectively excavating the locally underlying structure of the data that lies in the real world. Essentially, LFDA can be regarded as a weighted LDA by locally weighing and matrices. Therefore, the two modified scatter matrices, denoted as and , can be formulated as
(17)  
where the two weights ( and ) denote the samplewise similarities. There are several widelyused strategies in calculating such a similarity matrix symbolized by . A simple yet effective one is given by , if , where represents the nearestneighbor of ; otherwise, . Another commonlyused one was constructed based on the radial basis function (RBF) with a standard derivation of , as defined in (15). Please refer to [5, 36, 121] that might be useful for those who are interested in more types of .
Similar to SLDA that first projects the original data into a subspace and then LDA is performed in the transformed subspace, FSDA starts with maximizing the betweenspectral scatter matrix () to enhance the differences along the spectral dimension, and similarly, the LDA is further used for extracting the representations of class separability from the feature domain. In the first step, let be the average value of the th class and the th spectral band, then we have the definition of as follows:
(18) 
where is the spectral representation in the feature space and . The primary transformation () that aims at improving the spectral discriminant can be estimated with maximizing the trace term of as
(19) 
Using the obtained , the latent representation in the feature space can be further fed into the nextstep LDA.
GDA and Its Extensions: Before revisiting the GDA methods, we first introduce and formulate the general graph embedding (GGE) framework presented in [118] with Eq. (12). Obviously, the extracted features in the GGE framework are determined by the construction of to a great extent. Thus, we will highlight several kinds of representative affinity matrices corresponding to the different graph embedding approaches, i.e. LDA, LE [4] and its linearized LPP [35], LLE [95], sparse GDA (SGDA) [76], and collaborative GDA (CGDA) [75]. Fig. 8 visualizes the affinity matrices given by five different strategies in a fourclass case selected from Houston 2012 dataset.
LDAlike affinity matrix: In essence, LDA is vested in a special case of GGE framework with , whose affinity matrix can be represented as
(20) 
where is the number of samples belonging to th class.
LPP or LEbased affinity matrix: One is to be constructed in kernel space with a higher dimension via similarity measurement, i.e. extensively using (15).
LLEbased affinity matrix: Different from the handcrafted graph, LLE reconstructs each given sample with its nearest neighbors by exploiting the linear regression techniques [42, 37]. As a result, the reconstruction coefficients () can be obtained by solving the optimization problem of (10). With the known , it is straightforward to derive the needful affinity matrix, denoted as ,
(21) 
thereby inducing the Laplacian matrix as .
SGDA and CGDAguided affinity matrix: Similarly to LLE, the affinity matrix can be estimated using the datadriven representation learning, i.e., sparse and collaborative representations [125, 43, 38]. Accordingly, the two learning strategies can be equivalent to respectively solving the constrained norm optimization problem:
(22) 
and the norm optimization problem:
(23) 
The aforementioned affinity matrices can be unified to the GGE framework of (12).
In addition to SGDA and CGDA (the two baselines), Huang et al. [45] learned a set of sparse coefficients on manifolds and then preserved the sparse manifold structure in the embedded space. The work in [117] extended the existing SGDA to the spatialspectral graph embedding to address the issues of the spatial variability and spectral multimodality. With the embedding of the intrinsic geometric structure of the data, a Laplacian regularizer CGDA [69] was developed to further improve the graph’s confidence. Li et al. [70] simultaneously integrated the sparsity and lowrankness into the graph for capturing a more robust structure of the data locally and globally. Furthermore, Pan et al. [85] further improved the above work by unfolding the HSI data with the form of a tensor.
KDA: In reality, the HSI usually exhibits a highly nonlinear data distribution, which may result in difficulties to effectively identify the materials. The solution to this issue is making use of a socalled kernel trick [83] that can map the data of the input space into a new Hilbert space with a higher feature dimension. In the kernelinduced space, the complex nonlinearity of the HSI can be well analyzed in a linearized system. Comparatively, the input to KDA is an inner product of original data pairs, defined as which can be given by (15). By introducing the kernel Gram matrix with , most of previous LDAbased methods can be simply extended to the corresponding kernelized versions, i.e. KLDA and KLFDA can calculate their projections by solving a GED problem of
(24) 
Note that in KLDA, while and are computed by and in the kernel space, respectively, for KLFDA. Furthermore, for KSGDA and KCGDA, the main difference lies in the computation of the adjacency matrix, which can be performed in the kernel space by solving the general kernel coding problem as follows:
(25) 
where can be selected to be either sparsityprompting term of KSGDA or dense (or collaborative) term of KCGDA. In [22] and [69], the solutions in (25) have been theoretically guaranteed in the same way by solving the problems (22) and (23) using the alternating direction method of multiplier (ADMM) [9] and leastsquare regression with Tikhonov regularization [28], respectively.
Regressioninduced Representation Learning (RIRL)
RIRL provides a new insight from the regression’s point of view to model the FE behavior by bridging the training samples with the corresponding labels rather than indirectly using the label information in the form of graph or affinity matrix in DAFEbased methods.
LeastSquares Dimension Reduction (LSDR): We begin with sliced inverse regression (SIR) [66], which is a landmark in supervised FE techniques. It assumes that the pairwise data are conditionally independent on the tobeestimated subspace features , formulated as . Following this rule, the LSDR proposed by Suzuki and Sugiyama [106] attempts to find a maximizer of the squaredloss mutual information (SMI) to satisfy the previously mentioned independence assumption. The projections for LSDR can be searched by optimizing the following maximization problem:
(26) 
and the SMI to measure a statistical dependence between two discrete variables is defined as
(27) 
where is the probability distribution function.
LeastSquares Quadratic Mutual Information (LSQMI): Limited by the sensitivity of MI to outliers, authors of [96] designed a more robust LSQMI with the basis of QMI criterion, hence let us define the QMI as
(28) 
Similarly, we solve (26)like optimization problem by replacing SMI with QMI.
LeastSquares QMI Derivative (LSQMID): Due to the difficulty in accurately computing the derivative of QMI estimator, LSQMI was further extended to a computationally effective LSQMID by estimating the derivative of QMI instead of QMI itself [107]. In this work, authors have demonstrated a more accurate and efficient derivative computation of QMI.
Joint & Progressive Learning Strategy (JPlay): Another MIfree estimation group is latent subspace learning (LSL). One representative LSL performs FE and classification simultaneously in joint learning (JL) fashion [39]. With an expected output , the process can be modeled as
(29) 
where and are defined as the onehot encoded label matrices and the latent subspace projections, respectively. denotes the regression matrix that connects the learned subspace and the label information. can be formulated as
(30) 
In [39], the model’s solution has been proven to be a closedform. Moreover, the work in [40] explored a LDAlike graph as a regularizer to learn a spectrally discriminative feature representation, thus (29) becomes
(31)  
Beyond the JLbased models, Hong et al. [41] established a novel multilayered regression framework by following a joint and progressive learning strategy (JPlay). With the layerwise autoreconstruction mechanism effectively against spectral variabilities caused by complex noises and atmospheric effects, the linearized JPlay breaks through the performance bottleneck of traditional linear methods. More specifically, we have the resulting model in the following
(32)  
where the soft constraint can be used to relax the orthogonality. It is worth noting that such JLbased strategy can clearly tell the model which features are positive to classification task, owing to the joint strategy of feature extraction and classification.
Iv Deep Feature Extraction Techniques
Shallow feature extraction techniques often require careful engineering and domain knowledge of experts, which limit their applications. Different from them, deep learning techniques aim at automatically learning highlevel features from raw data in a hierarchical fashion. These features are more discriminative, abstract, and robust than shallow ones. Due to their powerful feature representation ability, deep learning techniques have been widely used to extract features from hyperspectral images in recent years [123, 67]. Among various deep learning models, autoencoders (AEs), convolutional neural networks (CNNs), and recurrent neural networks (RNNs), shown in Fig.9, are the most popular ones. In this section, we will present these models and their applications to hyperspectral feature extraction.
Iva AEs
As demonstrated in Fig.9, AE is mainly comprised of two modules: Encoder and decoder. Encoder maps the input vector into a hidden space , while decoder aims at getting a reconstruction result of the original input from . These processes can be formulated as
(33) 
where and denote the weights connecting the input layer to the hidden layer and the hidden layer to the output layer, respectively, and represent the biases of the hidden units and the output units, respectively, and is a nonlinear activation function. The training of AE is to minimize the residual between and . Once trained, the decoder is deleted and the hidden layer is considered as a feature representation of . In order to extract deep features, several AEs are often stacked together, generating a stacked AE (SAE) model. For SAE, the hidden layer in the preceding AE will be used as the input of the subsequent AE.
SAE is perhaps the earliest deep model used to extract features of hyperspectral images [72]. One typical benefit of SAEs is that each AE inside the network can be pretrained using both labeled and unlabeled samples, thus providing better initial values for network parameters compared to the random initialization. After layerwise pretraining, the finetuning of only a few layers can acquire satisfactory discriminant features. This training method is capable of alleviating overfitting problem when there only exist small numbers of training samples in hyperspectral images. In [14], the spectral information for each pixel was considered as a vector, and fed into an SAE model to extract deep features. These features extracted by SAEs can also be generalized from one image to another image, which was validated in [108] and [59].
In order to extract spatial features of each pixel, one often needs to select a local patch or cube centered at the pixel, and then input it into a feature extraction model. Since the inputs of SAEs are vectors, it is difficult to directly process patches or cubes. In [14] and [108], the local cubes from the first principal components of hyperspectral images were firstly reshaped into vectors, and then fed into SAEs to extract spatial features. In [56] and [16], Gabor features and extended morphological attribute profiles (i.e., the joint use of shallow and deep feature extraction methods) were used as the inputs of SAEs, making the network easier to extract highlevel spatial features. After the extraction of spectral and spatial features, they can be easily concatenated together to generate a spectralspatial joint feature [14, 108]. Compared to the concatenation method, the authors in [56] and [16] proposed to use another SAEs for fusing the spectral and spatial features, which may further enhance the discriminative ability of spectralspatial features.
Similar to the traditional feature extraction methods, one can also embed some prior or expected information into SAEs. Based upon the assumption that neighboring samples in the input space should have similar hidden representations, graph regularization was added to SAE for preserving such property [79, 105]. In [130], Zhou et al. imposed a local regularization via Fisher discriminant analysis on hidden layers to make the extracted features of samples from the same category close to each other while from different categories as far as possible, thus improving the discriminative ability of SAE. Meanwhile, they also added a diversity regularization term to make SAE extract compact features.
IvB Convolutional Neural Networks (CNNs)
CNNs are the most popularly adopted deep model for hyperspectral feature extraction. As shown in Fig.9, the basic components of a CNN model include convolutional layers, pooling layers, and fully connected layers. The convolutional layers are used to extract features with convolutional kernels (filters), which can be formulated as:
(34) 
where is the th feature maps, and denote the filters and biases of the th layer, respectively, and ’’ represents the convolutional operation. After the convolutional layer, the pooling layer is often adopted for reducing the size of the generated feature maps and producing more robust features. On the top of a CNN model, there often exist some fully connected layers, aiming at learning highlevel features and outputting the final results of the network.
For hyperspectral images, CNNs can be used for extracting spectral features [44] or spatial features [127, 13, 74], depending on the inputs of networks. In [44], Hu et al. designed a 1D CNN model to extract spectral features of each pixel. Compared to the traditional fullyconnected networks, CNNs have weightsharing and localconnection characteristics, making their training processes more efficient and effective. In [127], 2D CNN was explored to extract spatial features from a local cube. Different from SAEs, CNNs do not need to reshape the cube into a vector, thus preserving as much spatial information as possible. However, to make full use of the representation ability of CNNs, two important issues need to be considered. The first issue is the small number of training samples but highdimensional spectral information, which will easily lead to the overfitting problem. The second issue is the extraction of spectralspatial joint features, which can improve the classification performance in comparison with using the spectral or spatial feature only.
For the first issue, many widely used strategies in the field of natural image classification, such as dropout and weight decay, can be naturally adopted. In addition, a lot of promising methods have been proposed in the past few years. These methods can be divided into four different classes. The first class of methods is dimensionality reduction. In [127, 13, 102], PCA was employed to extract the first principal components of hyperspectral images as inputs of CNNs, thus simplifying the network structures. Similarly, a similaritybased band selection method was used in [78]. However, these dimensionality reduction methods are independent from the following CNNs, which may lose some useful information. Different from them, Ghamisi et al. proposed a novel method to adaptively select the most informative bands suitable for the CNN model [23]. The second class of methods is data augmentation. In [13], two methods were proposed to generate virtual samples. One is to multiply a random factor and add a random noise to training samples, while the other one is to combine two given samples from the same class with proper ratios. In [122], a data augmentation method based on distance density was proposed. Recently, Kong et al. proposed a random zero setting method to generate new samples [60]. The third class of methods is transfer learning. In [120] and [82], the authors found that CNNs trained by one hyperspectral data can be transferred to another data acquired by the same sensor, and finetuning only a few top layers achieves satisfying results. More interestingly, the works in [50, 15, 71] indicated that CNNs pretrained by natural images can be directly applied to extract spatial features of hyperspectral images. The fourth class of methods is semisupervised or even unsupervised learning. For examples, Wu and Prasad attempted to use a clustering model to obtain pseudo labels of unlabeled samples, and then combine the training samples and unlabeled samples (with their pseudo labels) together to train their network [113].
In terms of the second issue, one popularly used method is feeding a local cube, directly cropped from the original hyperspectral image, into a CNN with 3D convolution kernels for processing the spectral and spatial information simultaneously. The number of channels in the 3D convolutional kernel is smaller than or equal to that of its input layer. However, the former one dramatically increases the computational complexity due to the simultaneous convolution operators in both spectral domain and spatial domain, while the latter one heavily increases the number of parameters to optimize. Another candidate method is to decouple the task of spectralspatial feature extraction into two parts: Spectral feature extraction and spatial feature extraction. In [120] and [115], a parallel structure was employed to extract spectralspatial features. Specifically, 1D CNN and 2D CNN were designed to extract spectral features and spatial features,respectively; these two features were then concatenated together and fused via a few fullyconnected layers. Since 2D CNN focuses on extracting spatial features, some redundant spectral information can be preprocessed to reduce the computational complexity. In [128], a serial structure was also used to extract spectralspatial features. It firstly applied several convolutions to extract spectral features and then fed the extracted features into several 3D convolutions to extract spatial features.
IvC Recurrent Neural Networks (RNNs)
RNNs have been popularly employed to sequential data analysis, such as machine translation and speech recognition. Different from the feedforward neural network, RNN takes advantage of a recurrent edge for connecting the neuron to itself across time. Therefore, it is able to model the probability distribution of sequence data. To make this subsection easier to follow, we first provide a brief and general discussion on RNN. Then, we briefly describe how to use RNN specifically for the classification of hyperspectral images.
Fig.9 shows an example of RNN. Given a sequence , where generally denotes the information at the time , the output of the hidden layer at the th time step is:
(35) 
where and represent weight matrices from the current input layer to the hidden layer and the preceding hidden layer to the current hidden layer, respectively, is the output of the hidden layer at the preceding time, and is a bias vector. According to this equation, it can be observed that the contextual relationships in the time domain are constructed via a recurrent connection. Ideally, will capture most of the information, and can be considered as the final feature of the sequence data. In terms of classification tasks, one often inputs into an output layer , which can be described as:
(36) 
where is the weight matrix from the hidden layer to the output layer, and is a bias vector.
In recent years, RNNs have attracted more and more attention in the field of hyperspectral image feature extraction. To make full use of RNNs, one needs to first ask the following question: How to construct the sequence? An intuitive method is regarding the whole spectral bands as a sequence [73, 129]. For each pixel, its spectral values are fed into RNN from the first band to the last band, and the output of the hidden layer at the last band is the extracted spectral feature. Different from the traditional sequences in speech recognition or machine translation tasks, the succeeding bands do not depend on the preceding ones. Thus, Liu et al. also fed the spectral sequence from the last band to the first band to construct a bidirectional RNN model [73]. Another method is using a local patch or cube to construct the sequence [126, 100, 129]. For examples, Zhou et al. regarded the rows of each local patch, cropped from the first principal component of hyperspectral images, as a sequence, and fed them into RNN one by one to extract spatial features [129]; Zhang et al. adopted each pixel and its neighboring pixels in the cube to form a sequence [126]. These pixels were firstly sorted according to their similarities to the center pixel, and then fed into RNN sequentially to extract locally spatial features.
In real applications, the constructed sequence may be very long. Take the widely used Indian Pines data as an example, the length of the sequence is 200 (the number of spectral bands) if we use the first method mentioned above to construct the sequence. This sequence will increase the training difficulty, because of the gradients tending to either vanish or explode. In order to deal with this issue, long shortterm memory (LSTM) was employed as a more sophisticated recurrent unit [73, 77, 116, 129]. The core components of LSTM are three gates: an input gate, a forget gate, and an output gate. These gates together control the flow of information in the network. Similarly, gated recurrent unit (GRU) was also employed, which has only two gates (i.e., an update gate and a reset gate). Compared to LSTM units, GRUs have less numbers of parameters, which may be more suitable for hyperspectral image feature extraction since it usually has a limited number of training samples. Another candidate scheme to address the issue is to divide the longterm sequence into shorter sequences [116, 32]. For example, in [32], Hang et al. proposed to group the adjacent bands of hyperspectral images into subsequences, and then use RNNs to extract features from them. Since nonadjacent bands have some complementarity, they also used another RNN to fuse the extracted features.
IvD Integrated Networks
In general, AEs and RNNs are good at processing vectorized inputs, thus achieving promising results in terms of spectral feature extraction. However, both of them need to reshape the input patches or cubes into vectors during spatial feature extraction, which may destroy some spatial information. In contrast, CNNs are able to directly deal with image patches and cubes, resulting in more powerful spatial features than AEs and RNNs. It is natural to think whether we can integrate these networks together to make full use of their respective advantages. In the past few years, numerous works have been proposed in this direction.
One kind of integration method is to use each network independently, and then combining their results together [100, 116, 34, 32]. In [34], a parallel framework was proposed to extract spectralspatial joint features from hyperspectral images. In this framework, SAE was employed to extract spectral features of each pixel, and CNN was used to extract spatial features from the corresponding image patch. These two results were fused by a fully connected layer. Similar to this work, Xu et al. also adopted the parallel framework but used LSTM to extract the spectral features [116]. Different from them, Hang et al. proposed a serial framework to fuse CNNs and RNNs [32]. Specifically, they used CNN to extract the spatial features from each band of hyperspectral images, and then used RNN to fuse the extracted spatial features. In [100], Shi and Pun also employed a serial framework to integrate CNN and RNN for spectralspatial feature extraction.
Another kind of integration method is embedding the core component (i.e., convolutional operators) of CNNs into AEs or RNNs [59, 73]. In [59], an unsupervised spectralspatial feature extraction network was proposed. The whole framework was similar to AEs, which also adopted the socalled encoderdecoder paradigm. However, the fullyconnected operators in AEs were replaced by convolutional operators, so that the network can directly extract spectralspatial joint features from cubes. In [73], Liu et al. proposed a spectralspatial feature extraction method based on a convolutional LSTM network. Instead of fully connected operators, they also used convolutional operators in LSTM units. For a given cube, each band was fed into the convolutional LSTM unit sequentially. The convolutional operators can extract the spatial features while the recurrent operators can extract the spectral features. The whole network was optimized in an endtoend way, thus achieving satisfactory performance.
V Experimental Results
To evaluate the performance of different FE techniques, we selected four techniques from the UFE category (i.e., PCA [54], MSTV [17], OTVCA [90], LPP [35]), four techniques from the SFE category (i.e., LDA [18], CGDA [69], LSDR [106], and JPlay [41]), and five techniques from the deep FE category (i.e., SAE [111], RNN [32], CNN [61], CAE [58], and CRNN [73]). Here, we set the tuning parameters for those algorithms before representing the experimental results.
Va Algorithm Setup
The parameter setting usually plays a crucial role in assessing the performance of FE algorithms. Subspace dimension (or number of features, ) is a common parameter for all compared algorithms. Selection of the number of features is a hard task for hyperspectral image analysis. The endmember selection/extraction, subspace identification, and/or rank selection are all refer to this issue [27, 91, 8]. For a fair and simplified comparison, the parameter is unifiedly assigned to the same value, which is equal to the number of classes (). We should note that in LDA is automatically determined as , due to the class separability (fisher’s criterion).
Unsupervised FE

PCA: This method is a parameter free technique.

MSTV: In [17], all parameters are adjusted using a trial and error approach. The multiscale parameters adjusting the degree of smoothness (as suggested in [17]) are set to 0.003, 0.02, and 0.01. The spatial scale for the structure extraction in three levels (as suggested in [17]) is set to 2, 1, and 3.

OTVCA: This method is initialized as recommended in [90]. The tuning parameter which controls the level of smoothness applied on the features is set to one percent of the maximum intensity range of the datatsets.

LPP: The number of neighbors is set to 12. The bandwidth of the Gaussian kernel is set to 1.
Supervised FE
A common strategy for a model selection is to run crossvalidation (CV) on the training set, since the labeled samples are available in SFE. Therefore, we used the CV strategy on the following studied algorithms for parameter selection.

LDA: This method can be viewed as a baseline for SFE. There is no additional parameter in LDA.

CGDA: The Eq. (23) can be tuned to a regularized optimization problem, where one extra parameter – regularized norm – needs to be set in advance in the process of graph construction, which can be searched in the range of by CV. In the experiments, is used for all three datasets.

LSDR: Two parameters are involved in LSDR; The standard deviation for the Gaussian function and the regularization parameter which are selected in the range of and , respectively, using CV. Finally, and are both set to in our experiments.

JPlay: There are three regularization parameters (, , and ) that need to be set in the JPlay model (32). With the CV conducted on the training set of three different datasets, the regularization parameters are selected in the ranges of , yielding the final setting of for the first dataset, for the second dataset, and for the last dataset.
Deep FE

SAE: The input of SAE is the original spectral information of each pixel. Three hidden layers are used. The numbers of neurons from the first to the third hidden layer are set to 32, 64, and 128, respectively. ReLU is adopted as the activation function for each hidden layer.

CNN: The input of CNN is a small cube with a size of , where represents the number of spectral bands for each data. Three convolutional layers are used. Each convolutional layer is sequentially followed by a batch normalization layer, a ReLU activation function, and a maxpooling layer. Note that the last pooling layer is an adaptive maxpooling layer, making the output size equals to for any input sizes. The kernel size for each convolution is , and the numbers of kernels from the first to the third convolutions are set to 32, 64, and 128, respectively. Padding operators are used to preserve the spatial size after each convolutional operator.

PCNN: PCA is applied prior to CNN to reduce the spectral dimension of the HSI. The number of reduced dimension by PCA is set to the number of classes (). The input cube for CNN is of size .

RNN: The input of RNN is the same as the input of SAE. Two recurrent layers with GRU are employed. The number of neurons in each recurrent layer is set to 128.

Integrated Networks: Convolutional AE (CAE) and convolutional RNN (CRNN) are selected as two representative integrated networks. The input for them is the same as that for the CNN. For CAE, three convolutional layers and three deconvolutional layers are adopted. All of them use kernels. The numbers of kernels from the first to the third convolutional layers are set to 32, 64, and 128, respectively. In contrast, the numbers of kernels from the first and the third deconvolutional layers are set to 64, 32, and , respectively. Similar to [73], CRNN adopts two recurrent layers with convolutional LSTM units. For both recurrent layers, convolutional kernels are applied. The numbers of kernels for the first and the second recurrent layers are set to 32 and 64, respectively.
All of the above deep learning related models are implemented in the PyTorch framework. In order to optimize them, we use the Adam algorithm with default parameters. The batch size, the learning rate, and the number of training epochs are set to 128, 0.001, and 200, respectively. To reduce the effects of random initialization, all of the deep learning models are repeated 5 times, and the mean values are reported.
Random Forest Classifier
Apart from the deep FE technqiues, all the other FE techniques use random forest (RF) to perform the classification task. The number of trees selected for RF is set to 200. We approximately set the number of the prediction variable to the square root of the number of input bands.
Indian Pines 2010  
Shallow FE  Deep FE  
UFE  SFE  SFE  
Spectral  PCA  MSTV  OTVCA  LPP  LDA  CGDA  LSDR  JPlay  SAE  RNN  CNN  CAE  CRNN  PCNN  
1  0.9260  0.9064  0.9992  0.9260  0.8850  0.9628  0.8144  0.8459  0.9230  0.9327  0.8829  0.9275  0.9432  0.9590  0.9397 
2  0.8769  0.9976  1.0000  1.0000  0.9976  1.0000  0.9961  0.8933  0.9984  0.9573  0.9178  0.9544  0.9961  0.8596  0.9998 
3  0.8862  0.9724  0.9724  0.9897  0.9724  0.9724  0.9759  0.9655  0.9862  0.9683  0.9405  0.9890  0.9986  0.8241  0.9938 
4  0.6888  0.7762  1.0000  0.8953  0.8742  0.9270  0.7474  0.8194  0.8300  0.7508  0.7621  0.8840  0.8822  0.8569  0.8930 
5  0.8058  0.8855  0.8039  0.8151  0.8682  0.8802  0.8096  0.8394  0.8665  0.8488  0.8474  0.8692  0.8570  0.8638  0.8706 
6  0.8172  0.8797  0.9946  0.7094  0.9739  0.7883  0.8284  0.9397  0.9418  0.9013  0.9204  0.9268  0.9167  0.6975  0.9127 
7  0.4170  0.5954  0.6792  0.7166  0.6985  0.7166  0.6845  0.7030  0.6958  0.6265  0.5795  0.6507  0.6818  0.6348  0.7103 
8  0.2530  0.2583  0.2599  0.2768  0.2961  0.2955  0.2431  0.2952  0.2758  0.5349  0.2840  0.8776  0.6776  0.9934  0.4725 
9  0.6545  0.7498  0.7048  0.7943  0.8913  0.8452  0.7971  0.8419  0.8142  0.8732  0.8533  0.8302  0.8336  0.8194  0.8621 
10  0.8229  0.9406  0.9594  0.9368  0.9019  0.9804  0.8514  0.7761  0.9663  0.9015  0.8096  0.8368  0.8752  0.7946  0.9289 
11  0.6658  0.8402  0.9224  0.9945  0.8943  0.9288  0.7195  0.7651  0.8052  0.8121  0.7414  0.7440  0.7633  0.6492  0.9165 
12  0.9985  1.0000  1.0000  1.0000  0.9995  1.0000  1.0000  1.0000  1.0000  0.9998  0.9765  0.9945  0.9838  0.9748  0.9991 
13  0.9468  0.9962  0.9879  0.9888  0.9925  0.9830  0.9580  0.9738  0.9819  0.9621  0.9427  0.9925  0.9930  0.9892  0.9959 
14  0.8783  0.9000  0.9615  0.9145  0.9344  0.9174  0.8756  0.8628  0.8953  0.9094  0.9030  0.9984  0.9985  0.8981  0.9993 
15  0.9333  0.9667  0.9511  0.9933  0.9489  0.9756  0.9333  0.9311  0.9600  0.9307  0.8119  0.9947  0.9942  0.7556  0.9978 
16  0.3735  0.2036  0.2885  0.1601  0.5217  0.1719  0.3439  0.4901  0.4466  0.2053  0.2060  0.0980  0.1700  0.6028  0.1051 
AA  0.7465  0.8043  0.8428  0.8194  0.8532  0.8341  0.7861  0.8089  0.8367  0.8197  0.7737  0.8480  0.8478  0.8233  0.8498 
OA  0.7866  0.8598  0.8561  0.8378  0.9112  0.8748  0.8370  0.8748  0.8829  0.8836  0.8655  0.8945  0.8911  0.8525  0.9018 
0.7390  0.8297  0.8247  0.8054  0.8909  0.8481  0.8010  0.8466  0.8571  0.8580  0.8355  0.8716  0.8673  0.8213  0.8802  
Houston 2012  
Shallow FE  Deep FE  
UFE  SFE  SFE  
Spectral  PCA  MSTV  OTVCA  LPP  LDA  CGDA  LSDR  JPlay  SAE  RNN  CNN  CAE  CRNN  PCNN  
1  0.8262  0.8272  0.8025  0.8205  0.8110  0.8177  0.8139  0.8120  0.7768  0.8217  0.8182  0.8104  0.8154  0.8245  0.8089 
2  0.8318  0.8393  0.8412  0.8515  0.8214  0.8355  0.8327  0.8553  0.9662  0.8274  0.8153  0.8425  0.8167  0.8412  0.8293 
3  0.9782  1.0000  0.9822  1.0000  1.0000  1.0000  1.0000  1.0000  0.9980  0.9895  0.9939  0.8594  0.7731  0.9156  0.8432 
4  0.9138  0.9081  0.7633  0.8873  0.9479  0.8920  0.9053  0.8864  0.9564  0.9773  0.9040  0.9170  0.9153  0.9129  0.9159 
5  0.9659  0.9886  0.9915  0.9991  0.9867  0.9384  0.9915  0.9688  0.9782  0.9438  0.9389  0.9699  0.9585  0.9881  0.9824 
6  0.9930  0.9930  0.9580  0.9580  0.9790  1.0000  0.8741  0.9860  0.9930  0.9874  0.9678  0.8769  0.9776  0.9483  0.9497 
7  0.7463  0.8927  0.6362  0.7090  0.9123  0.7901  0.8535  0.8526  0.7817  0.7293  0.7392  0.8802  0.8694  0.8642  0.8627 
8  0.3305  0.4606  0.5992  0.6724  0.4311  0.7379  0.4302  0.4710  0.7806  0.3792  0.4153  0.6344  0.6762  0.5305  0.8351 
9  0.6771  0.7885  0.8706  0.9008  0.7413  0.6449  0.7186  0.6752  0.7592  0.7145  0.7367  0.8595  0.8540  0.8404  0.8691 
10  0.4295  0.4749  0.6612  0.8398  0.4595  0.4662  0.4826  0.5792  0.6014  0.5556  0.5373  0.5674  0.5782  0.4514  0.6168 
11  0.7011  0.7268  0.9820  0.9924  0.7306  0.7239  0.7287  0.5806  0.6983  0.6231  0.7250  0.7417  0.7292  0.6186  0.7913 
12  0.5485  0.9145  0.7349  0.9625  0.7560  0.6513  0.7656  0.5687  0.7858  0.6305  0.7606  0.9379  0.9402  0.8440  0.9593 
13  0.6140  0.7754  0.6982  0.7789  0.8105  0.6105  0.7719  0.6702  0.7509  0.4516  0.6656  0.8835  0.8968  0.8414  0.8765 
14  0.9838  0.9919  1.0000  1.0000  0.9960  0.9919  0.9879  0.9595  0.9879  0.9692  0.9850  0.9943  0.9773  0.9603  0.9968 
15  0.9789  0.9746  1.0000  0.9789  0.9746  0.9831  0.9852  0.9514  0.9831  0.9732  0.9607  0.8072  0.7471  0.9345  0.8592 
AA  0.7679  0.8371  0.8347  0.8901  0.8239  0.8056  0.8094  0.7878  0.8532  0.7716  0.7976  0.8388  0.8350  0.8210  0.8664 
OA  0.7278  0.8058  0.8088  0.8753  0.7874  0.7745  0.7789  0.7524  0.8280  0.7436  0.7646  0.8239  0.8184  0.7921  0.8526 
0.7076  0.7895  0.7923  0.8648  0.7700  0.7552  0.7604  0.7315  0.8134  0.7235  0.7469  0.8096  0.8036  0.7761  0.8404  
Houston 2017  
Shallow FE  Deep FE  
UFE  SFE  SFE  
Spectral  PCA  MSTV  OTVCA  LPP  LDA  CGDA  LSDR  JPlay  SAE  RNN  CNN  CAE  CRNN  PCNN  
1  0.3088  0.8781  0.0536  0.6842  0.6618  0.6256  0.7575  0.7969  0.5991  0.7940  0.5702  0.7516  0.4428  0.6338  0.6638 
2  0.7603  0.8396  0.7046  0.6376  0.8122  0.8474  0.8076  0.7747  0.8347  0.7893  0.6975  0.8173  0.8849  0.8707  0.8376 
3  1.0000  1.0000  1.0000  0.9972  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.9972  0.7739  0.8482  0.9924  0.8045 
4  0.9134  0.9494  0.6238  0.6775  0.9453  0.9059  0.9265  0.9276  0.9298  0.9221  0.8613  0.9444  0.9362  0.9439  0.9595 
5  0.4119  0.4668  0.2676  0.1679  0.4728  0.5258  0.4661  0.4289  0.3971  0.4982  0.4040  0.4330  0.5396  0.5404  0.4800 
6  0.2570  0.2990  0.3835  0.3164  0.3008  0.2910  0.2776  0.2726  0.2780  0.2585  0.2537  0.3050  0.3080  0.2902  0.3377 
7  0.3025  0.3025  0.3025  0.3025  0.3025  0.3025  0.3109  0.3025  0.2857  0.3025  0.3025  0.2908  0.2723  0.2997  0.3176 
8  0.7657  0.7675  0.7599  0.7417  0.7785  0.7849  0.7544  0.7518  0.7771  0.7216  0.7356  0.8538  0.8583  0.8092  0.8677 
9  0.3849  0.3877  0.5767  0.5990  0.4887  0.3917  0.3672  0.5255  0.5877  0.6302  0.4186  0.7970  0.7371  0.3717  0.8659 
10  0.3603  0.4360  0.3747  0.4491  0.4230  0.4086  0.3790  0.3497  0.4010  0.3819  0.3465  0.5902  0.4957  0.5484  0.5778 
11  0.4162  0.4792  0.7862  0.7596  0.5085  0.4667  0.4266  0.4422  0.5359  0.4143  0.4699  0.5456  0.5781  0.6048  0.5948 
12  0.0132  0.0046  0.0093  0.0077  0.0070  0.0023  0.0170  0.0000  0.0302  0.0152  0.0697  0.0511  0.0477  0.0927  0.0579 
13  0.4525  0.5556  0.4238  0.4090  0.5442  0.5164  0.5707  0.5603  0.5324  0.4523  0.4789  0.5148  0.5619  0.4246  0.5811 
14  0.3019  0.2629  0.5460  0.4665  0.3651  0.4152  0.2294  0.2073  0.3212  0.3789  0.3309  0.5289  0.6763  0.3375  0.5705 
15  0.6303  0.4721  0.4457  0.4887  0.4602  0.5549  0.4180  0.5234  0.5944  0.5197  0.5289  0.6277  0.6476  0.6447  0.6591 
16  0.6412  0.7611  0.6220  0.7648  0.7559  0.5888  0.7374  0.6403  0.6688  0.7457  0.7086  0.8498  0.7594  0.7173  0.8572 
17  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.9545  1.0000  1.0000  1.0000  1.0000  0.9545  0.8909  1.0000  1.0000 
18  0.4983  0.6885  0.6576  0.5197  0.7140  0.6581  0.6625  0.5686  0.7366  0.5346  0.6080  0.6365  0.5981  0.7692  0.7020 
19  0.5265  0.6363  0.9060  0.8777  0.7323  0.6989  0.6100  0.6266  0.6771  0.6569  0.6545  0.9102  0.8960  0.8006  0.9476 
20  0.4444  0.8904  0.9706  0.5253  0.8479  0.9189  0.6797  0.5955  0.5393  0.5388  0.4801  0.6246  0.5566  0.4879  0.6519 
AA  0.5195  0.6039  0.5707  0.5696  0.6060  0.5952  0.5676  0.5647  0.5863  0.5777  0.5458  0.6400  0.6268  0.6090  0.6667 
OA  0.4634  0.5101  0.5750  0.5899  0.5552  0.5027  0.4825  0.5492  0.5944  0.5938  0.4851  0.7278  0.6969  0.5116  0.7728 
0.3732  0.4317  0.4833  0.4974  0.4714  0.4231  0.4018  0.4560  0.5037  0.4948  0.3936  0.6474  0.6124  0.4372  0.7011  
(a) HSI  (b) LPP  (c) JPlay  (d) CNN  (e) PCNN  (f) Ground Reference 
(a) HSI  (b) OTVCA 
(c) JPlay  (d) CNN 
(e) PCNN  (f) Ground Reference 
(a) HSI  (b) OTVCA 
(c) JPlay  (d) CNN 
(e) PCNN  (f) Ground Reference 
VB Performance of the FE techniques on three HSIs
We have applied the FE techniques on the three hyperspectral datasets, i.e., Indian Pinse 2010, Houston 2012, and Houston 2017, and the classification accuracies including class accuracies, AA, OA, and coefficient are all shown in Tables V, VI, and VII, respectively. The results are first discussed within the categories and then between different categories. We should note that the results and discussions are in terms of classification accuracies obtained from the classification of the HSIs.
Unsupervised FE

PCA demonstrates the poorest performance compared with the other techniques, however, it considerably improves the classification accuracies compared with the results obtained by applying the RF on the spectral bands. One of the main disadvantages of the PCA is that, it does not take into account the noise and, therefore, the extracted features with having lower variance are often degraded by different types of noise existing in the hyperspectral image [93]. Additionally, PCA only takes into account the spectral correlation and it entirely neglects the spatial (neighboring) information.

LPP considerably outperforms the other UFE techniques for the Indian Pines dataset. However, in the case of Houston datasets, it provides very poor results. LPP incorporates the spatial information using the manifold learning process and by constructing the neighboring graph [35].

OTVCA outperforms the other UFE technqiues for Houston datasets. In the case of Houston 2012, the improvements are very considerable. OTVCA is robust to noise due to the signal model which takes into account the noise and model’s errors. Additionally, OTVCA exploits the spatial correlation by incorporating the TV penalty and, therefore, the extracted features are piecewise smooth and have high SNR [90]. Overall, it can be observed that OTVCA, which is a candidate from the lowrank reconstruction techniques, generally provides better classification accuracies than other the UFE techniques.
Supervised FE

LDA versus Spectral Classifier (RF): With the embedding of supervised information, LDA obviously performs better than the situation where RF is directly applied to the spectral signatures in terms of the overall performance and individual accuracies for most materials. This indicates the effectiveness of SFE to a great extent.

LDA versus CGDA: Although the classification performance of CGDA is inferior to that of LDA from the whole perspective, yet the advantage of CGDA mainly lies in its automation in computing the similarity (or connectivity) between samples. This could lead to a relatively stable FE, particularly in largescale and more complex hyperspectral scenes. Due to the datadriven graph embedding, CGDA yields a lower running speed than LDA in the process of model training.

LDA versus LSDR: Intuitively, LSDR provides competitive classification performance with LDA. However, LSDR is time consuming due to the distribution matching between input samples and labels. The requirement to estimate the statistical distribution also limits the LSDR’s stability, especially when the training set is available on a small scale (e.g., for the Indian Pines 2010 and Houston 2012 datasets).

LDA versus JPlay: Unlike the conventional regression techniques, JPlay is capable of extracting semantically meaningful and robust features, due to the multilayered structure and the selfreconstruction constraint (32). Quantitatively speaking, JPlay outperforms the other SFE methods. The CV provides a feasible solution to automatically determine the parameter combination in JPlay. Despite the ADMM solver designed for speeding up the optimization process such a multilayered parameter update inevitably suffers from high computational cost.
Deep FE

Spectral versus SpectralSpatial models: Most of the spectralspatial models (i.e., CNN, PCNN, CAE, and CRNN) achieve superior performance than spectral models (i.e., SAE and RNN) in terms of AA, OA, and Kappa due to the joint use of spectral and spatial information. This indicates that, besides the rich spectral information, spatial information is also important for hyperspectral image classification.

PCNN and CNN versus CAE and CRNN: Similar to SAE, CAE focuses on image reconstruction rather than classification. In contrast, PCNN and CNN are exclusively designed for the classification task, so they are able to learn more discriminative features than CAE, leading to better classification performance especially on the Houston 2017 dataset. Although CRNN also focuses on the classification task, it has a higher number of parameters to train. Using the same number of training samples and epochs, PCNN and CNN can achieve better results than CRNN in terms of AA, OA, and Kappa.

PCNN versus CNN: PCNN outperforms CNN in terms of the classification accuracies for all three datasets. We should note that the improvements are substantial in the case of the Houston 2012 and 2017 datasets. Due to the use of PCA, most of the redundant spectral information is reduced. Therefore, the number of trainable parameters in PCNN is smaller than that of CNN, making it easier to learn under the same condition.
Shallow UFE Versus Shallow SFE
For all three datasets used in the experiments, the UFE techniques provide better classification accuracies than the SFE techniques. Unlike the SFE, the UFE tends to pay more attentions on spatialspectral information extraction, owing to fully considering all samples of HSI as the model input. Conversely, the performance of SFE is, to a great extent, limited by the ability to largely gather HSI ground sampling. A direct evident is given in Tables V, VI, and VII. For the Indian Pines 2010 and Houston 2017 datasets where more training samples are available, SFEbased methods hold the competitive results with UFEbased ones, while for the Houston 2012 dataset, the classification performance of SFE is relatively inferior to that of UFE, due to the smallscale training set. Considering the low number of ground samples often available in HSI applications the experimental results confirm the advantage of UFE over SFE for HSI feature extraction.
Shallow FE versus Deep FE
At the first glance, the shallow FE approaches slightly outperform the deep FE techniques for the two datasets, i.e., Indian Pines 2010 and Houston 2012. However, a deep comparison reveals that some deep FE techniques, such as CNNbased FE, provide consistency and good performance over all three datasets. Additionally, when the dimensionreduced step (e.g., using PCA) is applied prior to the CNN technique, the resulting PCNN yields by far the second highest accuracies in the case of Indian Pines 2010 and Houston 2012 datasets (only moderately lower than LPP or OTVCA, respectively), and simultaneously obtains the best performance on the Houston 2017 dataset. It is worth mentioning that CNNbased FE methods obtain at least 10% increase over the shallow ones in the case of Houston 2017. This could be due to the high nonlinear behaviour of this dataset which contains 20 land cover classes. The main factors for CNNbased FE methods to obtain around 20% of improvement over those shallow FE methods on the Houston 2017 dataset are the availability of the sufficient training samples and modeling the spatial information of HSI well.
Comparisons of the Land Cover Maps
Figs. 10, 11, and 12 compare the classification maps for Indian Pines, Houston University 2012, and Houston University 2017, respectively. The figures compare the maps obtained from methods which provide the highest OA from each category (i.e., Shallow UFE, shallow SFE, and deep FE) along with the map obtained from the spectral classifier (HSI). Additionally, we depicted the maps obtained by CNN for all three datasets since it provides the highest OA among the deep FE techniques which do not exploit a reduction step.
Overall, the classification maps of either unsupervised or supervised FEbased approaches (e.g., LPP, JPlay, CNN, PCNN) are smoother compared to the results only using HSI that tends to generate sparse mislabeled pixels. More specifically, the classification maps generated by spectralspatial FEbased methods, e.g., OTVCA, CNN, and PCNN, are usually a bit oversmoothed, leading to the creation of fake structures, especially for the Indian Pines 2010 and Houston 2017 datasets. In the case of OTVCA, the over smoothing can be avoided by decreasing the tuning parameter. In contrast, JPlay obtains relatively desirable classification maps, despite the lack of spatial information modeling. It is worth mentioning that the JPlay algorithm can maintain the structural information for Houston 2012 in the shadow covered region where pixels at some bands are considerably attenuated. This is due to the elimination of the spectral variability using the selfreconstruction regularization (the third term in (32)) and the multilayered linearized regression technique.
VC Performance with respect to the Number of Training Samples
In this subsection, we investigate the performance of the FE techniques in terms of classification accuracies with respect to the number of training samples. As we have already stated, this analysis is of great interest due to two main reasons; first, ground sample acquisition and measurements are often cumbersome and could be impossible in cases for which the target area is not reachable. Additionally, the limited number of samples not only affects the performance of the supervised classifiers but also the supervised features extraction techniques since they are highly reliable on the number of training samples. Therefore, in this experiment, we perform an analysis on the Houston University dataset 2017 by comparing the performances of the feature extraction techniques by selecting 10, 25, 50, and 100 training samples randomly. Fig 13 compares the OAs obtained by applying RF on the spectral bands (labeled by HSI), the features extracted by OTVCA, and JPlay along with the OAs obtained by CNN, and PCNN. The results are mean values over 10 experiments by selecting the samples randomly (the standard deviations are shown by the error bars). The outcomes of the experiment can be summarized as follows:

The SFE technique (i.e., JPlay) improves the OAs compared to the spectral classifier. However, it provides much lower OA compared with the UFE and deep FE for all cases. Two aspects might explain this point. One is that JPlay fails to model spatial and contextual information, and another is that although JPlay attempts to enhance the reorientation ability of the features via multilayered linear mapping, yet it is still incomparable to the nonlinear deep FEbased techniques, particularly when the number of samples are increased.

In this experiment, the UFE technique (i.e. OTVCA) and the deep FE one, CNN, are performed similarly in terms of classification accuracies. Compared with the results given in Table VII, it can be observed that the random selection of the training samples over the entire class regions from the ground reference considerably improves the performance of RF applied on the features extracted by OTVCA. This is often due to the lack of a parameter selection technique to select the optimum parameter for the OTVCA algorithm which could lead to oversmoothing on the features.

The deep learning technique (i.e., PCNN) after using the reduction (i.e., PCA) provides very high OA for all the cases. Comparing the results with CNN (i.e., without using the PCA reduction) confirms the advantage of using the reduction stage prior to deep learning techniques.
Vi Conclusion and Summery
In the past decade, hyperspectral image FE has considerably evolved leading to three main research lines (i.e., shallow unsupervised, shallow supervised, and deep feature extraction approaches) that include the majority of feature extraction techniques, which were studied in this paper. The paper has systematically provided a technical overview of the stateoftheart techniques proposed in the literature by categorizing the aforementioned three focuses into subcategories. In order to make this research paper easytofollow for researchers at different levels (i.e., students, researchers, and senior researchers), we aimed at showing the evolution of each category over the decades rather than including many techniques with an exhaustive reference list. The experimental section was designed in a way to compare the performances of the techniques in two ways: (1) between all the categories (i.e., shallow unsupervised, shallow supervised, and deep feature extraction approaches) and (2) within each category by analyzing the corresponding subcategories. In this manner, a variety of subcategories is investigated detailing the evolution of shallow unsupervised approaches (i.e., conventional data projection schemes, band clustering/splitting techniques, lowrank reconstruction techniques, and manifold learning techniques), shallow supervised approaches (i.e., class separation discriminant analysis, graph embedding discriminant analysis, regressionbased representation learning, and joint & progressive learning strategy), and deep approaches (i.e., AE, CNN, RNN, and integrative approaches). Three recent hyperspectral datasets have been studied in this paper and the results are evaluated in terms of classification accuracies and the quality of the classification maps. The experiments carried out in this study showed that in terms of classification accuracies; 1) deep learning FE techniques (i.e., CNN and PCNN) can outperform the shallow ones in particular when a sufficient amount of training data is available 2) applying a dimensionality reduction step (such as PCA) prior to the deep learning technqiues considerably improves their performances 3) shallow UFE techniques not only outperform the SFE ones but also are very competitive compared with deep FE ones. However, we should mention that the conclusions are limited by the experiments carried out on the three HSI datasets. In addition, this paper provides an impressive amount of codes and libraries mostly written in Python and Matlab to ease out the task of researchers in this very vibrant field of research.
Acknowledgment
The authors would like to thank Prof. Melba Crawford for providing the Indian Pines 2010 Data and the National Center for Airborne Laser Mapping (NCALM), the University of Houston, and the IEEE GRSS Fusion Committee for providing the Houston datasets. The corresponding author of this paper is Dr. Danfeng Hong.
Footnotes
 https://www.emc.com/leadership/digitaluniverse/2014iview/executivesummary.htm
References
 (200303) Kernel independent component analysis. J. Mach. Learn. Res. 3, pp. 1–48. External Links: ISSN 15324435, Link, Document Cited by: §IIIA1.
 (200503) Exploiting manifold geometry in hyperspectral imagery. IEEE Trans. Geosci. Remote Sens. 43 (3), pp. 441–454. External Links: Document, ISSN 01962892 Cited by: §IIIA4.
 (2009) Classification of hyperspectral images with regularized linear discriminant analysis. IEEE Trans. Geos. Remote Sens. 47 (3), pp. 862–873. Cited by: §IIIB.
 (2002) Laplacian eigenmaps and spectral techniques for embedding and clustering. In Proc. NIPS, pp. 585–591. Cited by: §IIIA4, §IIIB1.
 (2003) Laplacian eigenmaps for dimensionality reduction and data representation. Neural Comput. 15 (6), pp. 1373–1396. Cited by: §IIIB1.
 (2015) Spectralspatial classification of hyperspectral remote sensing images. Artech House Publishers, INC, Boston, USA. Cited by: 1st item, §IA, §IC.
 (201204) Hyperspectral unmixing overview: geometrical, statistical, and sparse regressionbased approaches. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 5 (2), pp. 354–379. Cited by: §IIIA.
 (2008) Hyperspectral subspace identification. IEEE Trans. Geos. Remote Sens. 46 (8), pp. 2435–2445. Cited by: §VA.
 (201101) Distributed optimization and statistical learning via the alternating direction method of multipliers. Found. Trends Mach. Learn. 3 (1), pp. 1â122. External Links: ISSN 19358237, Link, Document Cited by: §IIIB1.
 (200210) Dimensionality reduction of hyperspectral data using discrete wavelet transform feature extraction. IEEE Trans. Geosci. Remote Sens. 40 (10), pp. 2331–2338. External Links: Document, ISSN 01962892 Cited by: §IIIA1.
 (2016) Supervised band selection using local spatial information for hyperspectral image. IEEE Geosci. Remote Sens. Lett. 13 (3), pp. 329–333. Cited by: §IIIB.
 (201105) BandClust: an unsupervised band reduction method for hyperspectral remote sensing. IEEE Geosci. Remote Sens. Lett. 8 (3), pp. 565–569. External Links: Document, ISSN 1545598X Cited by: §IIIA2.
 (2016) Deep feature extraction and classification of hyperspectral images based on convolutional neural networks. IEEE Trans. Geosci. Remote Sens. 54 (10), pp. 6232–6251. Cited by: §IVB, §IVB.
 (2014) Deep learningbased classification of hyperspectral data. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 7 (6), pp. 2094–2107. Cited by: §IVA, §IVA.
 (201811) Exploring hierarchical convolutional features for hyperspectral image classification. IEEE Trans. Geosci. Remote Sens. 56 (11), pp. 6712–6722. Cited by: §IVB.
 (2019) Active transfer learning network: a unified deep joint spectral–spatial feature learning model for hyperspectral image classification. IEEE Trans. Geosci. Remote Sens. 57 (3), pp. 1741–1754. Cited by: §IVA.
 (2019) Noiserobust hyperspectral image classification via multiscale total variation. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 12 (6), pp. 1948–1962. Cited by: §IIIA2, 2nd item, §V.
 (1973) Pattern classification and scene analysis. Vol. 3, Wiley New York. Cited by: §IIIB, §V.
 (2009) Kernel principal component analysis for the classfication of hyperspectral remotesensing data over urban areas. EURASIP Jour. Adv. Signal Proc., pp. 1–14. Cited by: §IIIA1.
 (2013) Advances in spectralspatial classification of hyperspectral images. Proceedings of the IEEE 101 (3), pp. 652–675. Cited by: §IC.
 (1990) Introduction to statistical pattern recognition.. 2nd edition, Academic Press, Inc., San Diego, CA. Cited by: §IA.
 (2010) Kernel sparse representation for image classification and face recognition. In Proc. ECCV, pp. 1–14. Cited by: §IIIB1.
 (201610) A selfimproving convolution neural network for the classification of hyperspectral data. IEEE Geosci. Remote Sens. Lett. 13 (10), pp. 1537–1541. External Links: ISSN 1545598X Cited by: §IVB.
 (201505) A survey on spectralâspatial classification techniques based on attribute profiles. IEEE Trans. Geosci. Remote Sens. 53 (5), pp. 2335–2353. Cited by: §IIIA.
 (2018Sep.) New frontiers in spectralspatial hyperspectral image classification: the latest advances based on mathematical morphology, markov random fields, segmentation, sparse representation, and deep learning. IEEE Geosci. Remote Sens. Mag. 6 (3), pp. 10–43. External Links: Document, ISSN 23737468 Cited by: §IC.
 (201703) Advanced spectral classifiers for hyperspectral images: a review. IEEE Geosci. Remote Sens. Mag. 5 (1), pp. 8–32. External Links: Document, ISSN 23737468 Cited by: §I.
 (201712) Advances in hyperspectral image and signal processing: a comprehensive overview of the state of the art. IEEE Geos. Remote Sens. Mag. 5 (4), pp. 37–78. Cited by: §VA.
 (1999) Tikhonov regularization and total least squares. SIAM Journal on Matrix Analysis and Applications 21 (1), pp. 185–194. External Links: Document, Link, https://doi.org/10.1137/S0895479897326432 Cited by: §IIIB1.
 (1988) A transformation for ordering multispectral data in terms of image quality with implications for noise removal. IEEE Trans. Geos. Remote Sens. 26 (1), pp. 65–74. Cited by: §IIIA1.
 (201507) The EnMAP spaceborne imaging spectroscopy mission for earth observation. Remote Sens. 7 (7), pp. 8830–8857. Cited by: §IB.
 (2006) Band selection for hyperspectral image classification using mutual information. IEEE Geosci. Remote Sens. Let. 3 (4), pp. 522–526. Cited by: §IIIB.
 (2019) Cascaded recurrent neural networks for hyperspectral image classification. IEEE Trans. Geosci. Remote Sens. 57 (8), pp. 5384–5394. Cited by: §IVC, §IVD, §V.
 (2016) Matrixbased discriminant subspace ensemble for hyperspectral image spatial–spectral feature fusion. IEEE Trans. Geos. Remote Sens. 54 (2), pp. 783–794. Cited by: §IIIB.
 (2018) Twostream deep architecture for hyperspectral image classification. IEEE Trans. Geosci. Remote Sens. 56 (4), pp. 2349–2361. Cited by: §IVD.
 (2003) Locality Preserving Projections. In Proc. NIPS, S. Thrun, L. Saul and B. Scholkopf (Eds.), Cited by: §IIIA4, §IIIB1, item , §V.
 (2004) Locality preserving projections. In Proc. NIPS, pp. 153–160. Cited by: §IIIB1.
 (2019) Learning to propagate labels on graphs: an iterative multitask regression framework for semisupervised hyperspectral dimensionality reduction. ISPRS J. Photogramm. Remote. Sens. 158, pp. 35–49. Cited by: §IIIB1.
 (2019) An augmented linear mixing model to address spectral variability for hyperspectral unmixing. IEEE Trans. Image Process. 28 (4), pp. 1923–1938. Cited by: §IIIB1.
 (2019) Cospace: common subspace learning from hyperspectralmultispectral correspondences. IEEE Trans. Geos. Remote Sens. 57 (7), pp. 4349–4359. Cited by: §IIIB2.
 (2019) Learnable manifold alignment (LeMA): a semisupervised crossmodality learning framework for land cover and land use classification. ISPRS J. Photogramm. Remote. Sens. 147, pp. 193–205. Cited by: §IIIB2.
 (2018) Joint & progressive learning from highdimensional data for multilabel classification. In Proc. ECCV, pp. 469–484. Cited by: §IIIB2, §V.
 (2017) Learning a robust local manifold representation for hyperspectral dimensionality reduction. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 10 (6), pp. 2960–2975. Cited by: §IIIB1.
 (2018) SULoRA: subspace unmixing with lowrank attribute embedding for hyperspectral data analysis. IEEE J. Sel. Topics Signal Process. 12 (6), pp. 1351–1363. Cited by: §IIIB1.
 (2015) Deep convolutional neural networks for hyperspectral image classification. Sensors 2015. Cited by: §IVB.
 (2015) Dimensionality reduction of hyperspectral images based on sparse discriminant manifold embedding. ISPRS J. Photogramm. Remote. Sens. 106, pp. 42–54. Cited by: §IIIB1.
 (1968) On the mean accuracy of statistical pattern recognizers. IEEE Trans. Inf. Theory 14, pp. 55 – 63. Cited by: §IA.
 (2001) Independent component analysis. Adaptive and Learning Systems for Signal Processing, Communications and Control Series, Wiley. External Links: ISBN 9780471405405, LCCN 2001024709 Cited by: §IIIA1.
 (2015) Feature space discriminant analysis for hyperspectral data feature reduction. ISPRS J. Photogramm. Remote. Sens. 102, pp. 1–13. Cited by: §IIIB.
 (201303) Feature mining for hyperspectral image classification. Proc. IEEE 101 (3), pp. 676–697. Cited by: §ID.
 (201710) Deep fully convolutional networkbased spatial distribution prediction for hyperspectral image classification. IEEE Trans. Geosci. Remote Sens. 55 (10), pp. 5585–5599. Cited by: §IVB.
 (1998) Supervised classification in highdimensional space: geometrical, statistical, and asymptotical properties of multivariate data. IEEE Trans. Syst. Man Cybern. C: Appl. Rev. 28 (1), pp. 39–54. Cited by: §IA.
 (2019Jan.) Intrinsic image recovery from remote sensing hyperspectral images. IEEE Trans. Geosci. Remote Sens. 57 (1), pp. 224–238. External Links: Document, ISSN 01962892 Cited by: §IIIA2.
 (2017Aug.) Superpixelbased intrinsic image decomposition of hyperspectral images. IEEE Trans. Geosci. Remote Sens. 55 (8), pp. 4285–4295. External Links: Document, ISSN 01962892 Cited by: §IIIA2.
 (1986) Principal component analysis. Springer Verlag. Cited by: §IIIA1, §V.
 (2014June) Feature extraction of hyperspectral images with image fusion and recursive filtering. IEEE Trans. Geos. Remote Sens. 52 (6), pp. 3742–3752. External Links: Document, ISSN 01962892 Cited by: §IIIA2.
 (2018) Classification of hyperspectral images by gabor filtering based deep network. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 11 (4), pp. 1166–1178. Cited by: §IVA.
 (2015April) Intrinsic image decomposition for feature extraction of hyperspectral images. IEEE Trans. Geosci. Remote Sens. 53 (4), pp. 2241–2253. External Links: Document, ISSN 01962892 Cited by: §IIIA2.
 (2017) Selftaught feature learning for hyperspectral image classification. IEEE Trans. Geosci. Remote Sens. 55 (5), pp. 2693–2705. Cited by: §V.
 (2017) Selftaught feature learning for hyperspectral image classification. IEEE Trans. Geosci. Remote Sens. 55 (5), pp. 2693–2705. Cited by: §IVA, §IVD.
 (201811) Spectralspatial feature extraction for hsi classification based on supervised hypergraph and sample expanded cnn. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 11 (11), pp. 4128–4140. Cited by: §IVB.
 (2012) Imagenet classification with deep convolutional neural networks. In Proc. NIPS, Lake Tahoe, Nevada, USA, pp. 1527–1554. Cited by: §V.
 (200107) Bestbases feature extraction algorithms for classification of hyperspectral data. IEEE Trans. Geosci. Remote Sens. 39 (7), pp. 1368–1379. External Links: Document, ISSN 01962892 Cited by: §IIIA2.
 (2003) Signal theory methods in multispectral remote sensing. edition, John Wiley & Sons, Hoboken, NJ. External Links: ISBN Cited by: §IA.
 (1993) Feature extraction based on decision boundaries. IEEE Trans. Pat. Analysis Mach. Intell. 15 (4), pp. 388–400. Cited by: §IIIB.
 (1990) Enhancement of high spectral resolution remotesensing data by a noiseadjusted principal components transform. IEEE Trans. Geos. Remote Sens. 28 (3), pp. 295–304. Cited by: §IIIA1.
 (1991) Sliced inverse regression for dimension reduction. J. Am. Stat. Assoc. 86 (414), pp. 316–327. Cited by: §IIIB2.
 (2019Sep.) Deep learning for hyperspectral image classification: an overview. IEEE Trans. Geosci. Remote Sens. 57 (9), pp. 6690–6709. External Links: Document, ISSN 15580644 Cited by: §IV.
 (201803) Discriminant analysisbased dimension reduction for hyperspectral image classification: a survey of the most recent advances and an experimental comparison of different techniques. IEEE Geosci. Remote Sens. Mag. 6 (1), pp. 15–34. External Links: Document, ISSN 23737468 Cited by: §ID.
 (2016) Laplacian regularized collaborative graph for discriminant analysis of hyperspectral imagery. IEEE Trans. Geosci. Remote Sens. 54 (12), pp. 7066–7076. Cited by: §IIIB1, §IIIB1, §V.
 (2016) Sparse and lowrank graph for discriminant analysis of hyperspectral imagery. IEEE Trans. Geosci. Remote Sens. 54 (7), pp. 4094–4105. Cited by: §IIIB1.
 (201808) Deep multiscale spectralspatial feature fusion for hyperspectral images classification. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 11 (8), pp. 2911–2924. Cited by: §IVB.
 (2018) Spectral transformation based on nonlinear principal component analysis for dimensionality reduction of hyperspectral images. Eur. J. Remote Sens. 51 (1), pp. 375–390. Cited by: §IVA.
 (2017) Bidirectionalconvolutional lstm based spectralspatial feature learning for hyperspectral image classification. Remote Sens. 9 (12), pp. 1330. Cited by: §IVC, §IVC, §IVD, 5th item, §V.
 (2019) Stfnet: a twostream convolutional neural network for spatiotemporal image fusion. IEEE Trans. Geosci. Remote Sens. 57 (9), pp. 6552–6564. Cited by: §IVB.
 (2014) Collaborative graphbased discriminant analysis for hyperspectral imagery. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 7 (6), pp. 2688–2696. Cited by: §IIIB1.
 (2014) Sparse graphbased discriminant analysis for hyperspectral imagery. IEEE Trans. Geosci. Remote Sens. 52 (7), pp. 3872–3884. Cited by: §IIIB1.
 (2019) Hyperspectral image classification using similarity measurementsbased deep recurrent neural networks. Remote Sens. 11 (2), pp. 194. Cited by: §IVC.
 (2018) Hyperspectral image classification based on deep deconvolution network with skip architecture. IEEE Trans. Geosci. Remote Sens. 56 (8), pp. 4781–4791. Cited by: §IVB.
 (2016) Spectral–spatial classification of hyperspectral image based on deep autoencoder. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 9 (9), pp. 4073–4085. Cited by: §IVA.
 (200907) ICA and kernel ICA for change detection in multispectral remote sensing images. In Proc. IGARSS, Vol. 2, pp. II–980–II–983. External Links: Document, ISSN 21536996 Cited by: §IIIA1.
 (200712) Clusteringbased hyperspectral band selection using information measures. IEEE Trans. Geosci. Remote Sens. 45 (12), pp. 4158–4171. External Links: Document, ISSN 01962892 Cited by: §IIIA2.
 (201708) Learning sensorspecific spatialspectral features of hyperspectral images via convolutional neural networks. IEEE Trans. Geosci. Remote Sens. 55 (8), pp. 4520–4533. Cited by: §IVB.
 (2001) An introduction to kernelbased learning algorithms. IEEE Trans. Neural Netw. 12 (2). Cited by: §IIIB1.
 (2011Mar) Kernel maximum autocorrelation factor and minimum noise fraction transformations. IEEE Trans. Image Process. 20 (3), pp. 612–624. Cited by: §IIIA1.
 (2017) Hyperspectral dimensionality reduction by tensor sparse and lowrank graphbased discriminant analysis. Remote Sens. 9 (5), pp. 452. Cited by: §IIIB1.
 (2015)