Manifold regularization for sparse unmixing of hyperspectral images

Background Recently, sparse unmixing has been successfully applied to spectral mixture analysis of remotely sensed hyperspectral images. Based on the assumption that the observed image signatures can be expressed in the form of linear combinations of a number of pure spectral signatures known in advance, unmixing of each mixed pixel in the scene is to find an optimal subset of signatures in a very large spectral library, which is cast into the framework of sparse regression. However, traditional sparse regression models, such as collaborative sparse regression, ignore the intrinsic geometric structure in the hyperspectral data. Results In this paper, we propose a novel model, called manifold regularized collaborative sparse regression, by introducing a manifold regularization to the collaborative sparse regression model. The manifold regularization utilizes a graph Laplacian to incorporate the locally geometrical structure of the hyperspectral data. An algorithm based on alternating direction method of multipliers has been developed for the manifold regularized collaborative sparse regression model. Conclusions Experimental results on both the simulated and real hyperspectral data sets have demonstrated the effectiveness of our proposed model.

endmembers, and their corresponding fractional abundances, has been proposed and widely applied in hyperspectral remote sensing data analysis (Bioucas-Dias and Plaza 2013).
For spectral unmixing, two kinds of models-the linear mixture model (LMM) (Singer et al. 1979) and the non-linear mixture model (NMM) (Hapke 1981)-have been proposed to characterize the mixed pixels. Although the LMM is not as accurate as the NMM to capture the mixing behavior of mixed pixels, it is more popular than the NMM for solving the spectral unmixing problem because of its simplicity and efficiency in most cases (Fan et al. 2009). In addition, the rapidly developed methods in classical signal processing field also provide effective tools to the solution of the LMM (Ma et al. 2014). Nevertheless, modeling the mixed pixels is a very complex and difficult task. In practice, we have to make a compromise between model accuracy and tractability. Therefore, we here focus on the LMM in this study.
The standard LMM used for spectral unmixing assumes that each pixel spectrum is a linear combination of the endmembers present in the scene weighted by the corresponding abundances. From the convex geometry point of view, the LMM forces the mixed pixels to belong to a simplex (or a convex hull), and the vertices of simplex correspond to the endmembers. Based on the geometrical interpretation, many spectral unmixing algorithms have been proposed for endmember extraction such as the N-FINDR (Winter 1999), pixel purity index (PPI) (Boardman et al. 1995), vertex component analysis (VCA) (Nascimento and Bioucas-Dias 2005), simplex growing algorithm (SGA) (Chang 2006) and their variants (Chan et al. 2011;Liu and Zhang 2012;Chang et al. 2010), and for abundance estimation such as the fully constrained least squares (FCLS) (Heinz and Chang 2001), distance geometry-based abundance estimation (DGAE) (Pu et al. 2014), and so on. However, these geometrical-based algorithms are likely to fail when the pixels are highly mixed. As an alternative, the statistical algorithms have been developed by formulating the spectral unmixing as a statistical inference problem. Such algorithms include the dependent component analysis (DECA) (Nascimento and Bioucas-Dias 2012), beta compositional model (BCM) (Zare et al. 2013) and normal compositional model (NCM) (Stein 2003) methods. Although the statistical algorithms have a natural framework for incorporating various priors and endmember variability (Somers et al. 2011;Zare and Ho 2014), it is difficult to derive the close-form expressions of the inference parameters and thus they suffer from high computational complexity.
Most of the spectral unmixing algorithms based on the standard LMM can not automatically determine the number of endmembers present in the scene. In addition, some endmembers produced by these algorithms are not necessarily present in the image, producing the so-called virtual endmembers (Chen 2011). The virtual endmembers can compensate the approximation of the LMM but will result in unidentifiability of the materials. To tackle these problem, the standard LMM has been extended into a semisupervised version (Liu and Zhang 2014;Iordache et al. 2012;Zhong and Zhang 2014;Feng et al. 2014;Iordache et al. 2011), i.e., by assuming that the endmembers are known in advance. Typically, Iordache et al. (2011) have proposed the sparse regression (SR) model by assuming that the endmembers present in the scene belong to a subset of samples available a priori in a library. The unmixing based on SR is called sparse unmixing. Experimental results have illustrated the potential of sparse unmixing in abundance estimation. The success of sparse unmixing relies crucially on the availability of suitable hyperspectral libraries because the libraries are hardly acquired under the same conditions of the remotely sensed images. Fortunately, this problem can be overcome by a delicate calibration procedure to adapt the library to the image or learning of the libraries directly from the data set without other priori information (Charles et al. 2011). The SR problem can be efficiently solved via the sparse unmixing algorithm via variable splitting and augmented Lagrangian (SUnSAL) (Iordache et al. 2011;Bioucas-Dias and Figueiredo 2010) by exploiting the sparse prior induced by the ℓ 1 norm. However, a high correlation of the spectral signatures limits the unmixing accuracy. To mitigate this limitation, Iordache et al. (2014a) have developed a collaborative SR (CSR) model by considering the structured sparsity, which exploits the fact that only a few spectral signatures in the library are active, in other words, only a few lines of abundances collected in a matrix are nonzero. Some modifications of the CSR can be found in Iordache et al. (2014b), Tang et al. (2015). However, the improvements in Iordache et al. (2014b), Tang et al. (2015) are limited since only the spectral information is considered to estimated the abundances.
Generally, the size of spectral library is often large, while the number of endmembers present in the scene is very small. Therefore, the fractional abundances are more likely to reside on a low-dimensional submanifold of the high-dimensional ambient Euclidean space. However, existing sparse unmixing methods only consider the Euclidean structure of the data space while ignoring the intrinsic manifold structure of the hyperspectral data. Many previous studies (Lu et al. 2013;Zheng et al. 2011;Guan et al. 2011;Seung and Lee 2000;He and Niyogi 2004;Belkin et al. 2006) have shown that exploiting the local geometrical structure (i.e., the intrinsic manifold structure) is very important to the model learning and data representation. In this paper, we incorporate the manifold regularization (Belkin et al. 2006) to the CSR and develop a novel model, called manifold regularized collaborative sparse regression (MCSR) model. The manifold regularization is characterized by a Laplacian graph which captures the local geometrical structure of the data manifold such that nearby mixed pixels in the intrinsic geometry of the data space are likely to have similar fractional abundances. By adding an additional manifold structure learning term to CSR, our proposed MCSR model is expected to have higher unmixing accuracy than CSR. To solve the MCSR, an optimization algorithm based on alternating direction method of multipliers (ADMM) is developed. It should be noted that similar works of using the manifold regularization have also been introduced in Lu et al. (2013), Tong et al. (2014), but they are different with ours in the following two main aspects. First, the works in Lu et al. (2013), Tong et al. (2014) are based on the standard LMM while our model is an extension of the SR model. Second, multiplicative iterative algorithm is used to optimize the nonnegative matrix factorization model in Lu et al. (2013), Tong et al. (2014), whereas the ADMM is used to cope with the proposed MCSR model.

Related work
In this section, we first describe the sparse unmixing problem and then briefly review the CSR model.

Sparse unmixing
The standard LMM signors the intimate mixture and holds at a macroscopic level (Bioucas-Dias and Plaza 2012). For a pixel spectrum y ∈ R l with l spectral bands, LMM assume that it is a linear combination of the endmembers present in the scene weighted by the corresponding abundances, i.e., where α = (α 1 , . . . , α p ) T ∈ R p is the abundance vector, M ∈ R l×p is the endmember matrix with each column m i ∈ R l being an endmember present in the scene, e ∈ R l denotes the noise or error term, and p is the number of endmembers present in the scene. To be physically meaningful, the abundance vector is usually subject to (s.t.) the sum-to-one and nonnegativity constraints Over the past decades, lots of spectral unmixing algorithms (Keshava 2003;Plaza 2012, 2013) have been developed based on the above LMM. However, the abundance estimation by these algorithms usually relies on the availability of pure spectral signatures in the input data or on their capacities of extracting endmembers. In addition, some algorithms perform unmixing by assuming all of the endmembers in M are present in the scene and by exploiting the sparsity prior of abundance α. If the abundance α is not sparse or sufficiently sparse, the results obtained by these algorithms will not be as accurate as we expect. To cope with these problems, sparse unmixing (Iordache et al. 2011) has been introduced based on the assumption that the endmember set {m 1 , m 2 , . . . , m p } present in the scene is contained in a spectral library denoted by {a 1 , a 2 , . . . , a m } known in advance, i.e., {m 1 , m 2 , . . . , m p } ⊂ {a 1 , a 2 , . . . , a m }. With the ever-growing availability of spectral libraries, the number of endmembers present in the scene is much less than the total number of endmembers in spectral library; thus, we have p ≪ m. By this way, unmixing is to find the optimal subset of signatures for the mixed pixels in the spectral library. The sparse unmixing model can be written as where A = [a 1 , a 2 , . . . , a m ] ∈ R l×m is the spectral library, and x ∈ R m is the abundance vector corresponding the spectral library A. Clearly, x is sparse. By exploiting the sparse prior of abundance x through the well-known ℓ 1 norm, we can estimate abundance vector x by the following sparse regression (SR) problem: where �x� 1 = m i=1 |x i | and SR > 0 is a sparse regularization parameter. It is worth to note that the sum-to-one constraint is no longer necessary since the nonnegativity constraint can automatically lead to the sum-to-one constraint as stated in Bioucas-Dias and Plaza (2012), Iordache et al. (2011). However, the high mutual coherence of the spectral library A limits the unmixing accuracy.

Collaborative sparse regression
In fact, sparse unmixing, as a semi-supervised model, is a typical underdetermined linear system. To solve it, sparsity prior for the fractional abundance of each individual pixel is imposed in (6). Although the high level of sparsity of fractional abundances can enhance the recovery ability of the ℓ 1 -minimization problem (6), the highly correlated samples in library still restrict the capability of model (6) to obtain a desirable solution.
To tackle this problem, the CSR model has been recently proposed in Iordache et al. (2014a). Unlike the ℓ 1 -minimization problem (6), the CSR model simultaneously (or collaboratively) imposes a sparsity to all pixels in the data set by exploiting the fact that pixels in a scene should share the same set of active endmembers, and thus only a few rows of the abundance matrix are nonzero.
Let Y = [y 1 , y 2 , . . . , y n ] ∈ R l×n be the n observed pixels arranged in a matrix, we can rewrite the sparse unmixing model (5) in the matrix form as where E = [e 1 , e 2 , . . . , e n ] ∈ R l×n is the noise matrix and X = [x 1 , x 2 , . . . , x n ] ∈ R m×n is the fractional abundance matrix. With the collaborative sparsity induced by the rowsparsity regularizer ℓ 2,1 norm, the CSR problem can be formulated as where �X� 2,1 = m i=1 �x i � 2 , x i represents the ith row of abundance fraction matrix X and CSR is a regularization parameter. An algorithm called collaborative sparse unmixing via variable splitting and augmented Lagrangian (CLSUnSAL) is provided in Iordache et al. (2014a) and a number of experiments have shown that the imposed collaborative sparsity prior can significantly reduce the probability of recovery failure.

Methods
In this section, we introduce an enhanced CSR model, called manifold regularized CSR (MCSR) model, by incorporating a manifold regularization to CSR, and then an alternating direction method is developed to solve the resulting optimization problem.

MCSR
As for CSR problem (8), the data fitting term �Y − AX� 2 F is useful for learning the Euclidean structures in the hyperspectral data space. As we have previously mentioned, the size of spectral library is usually very large, while the number of endmembers present in the scene is very small. From a geometric viewpoint, the fractional abundances are more likely to reside on a low-dimensional submanifold of the high-dimensional ambient Euclidean space. Recent studies (Lu et al. 2013;Zheng et al. 2011;Guan et al. 2011;Seung and Lee 2000; He and Niyogi 2004;Belkin et al. 2006) have shown that intrinsic geometric structures on manifolds are very important to the data representation and many manifold learning methods (Tenenbaum et al. 2000;Roweis and Saul 2000;Donoho and Grimes 2003;Belkin and Niyogi 2003;Lin and Zha 2008) have been proposed to recover the geometry of a data set. In the literature of spectral unmixing, many existing methods (Lu et al. 2013) only explore the Euclidean structure while fail to discover the intrinsic geometry structure of the data manifold. Therefore, we want to explore the ability of the intrinsic geometry structure of the hyperspectral data in improving the unmixing accuracy.
To preserve the intrinsic geometry structure of the data, a natural assumption, referred to the manifold assumption (He and Niyogi 2004), is that nearby data points are also nearby points in their low-dimensional representations. However, modeling the global geometric structures of the data is a very big challenge due to the insufficient number of samples and the high dimensionality of the ambient space. In practice, a nearest neighbor graph on the data points is often used to characterize the underlying local geometric structures.
Given n data points {y 1 , y 2 , . . . , y n } ∈ R l sampled from the underlying submanifold, we construct a nearest neighbor graph G with its ith node corresponding to the data point y i , i = 1, 2, . . . , n. For each node y i , one can put an edge between it and its k nearest neighbors. Let N k (y i ) = {y 1 i , . . . , y k i } be the set of its k nearest neighbors. And we define the weight matrix W on the graph G as The weight w ij is used to measure the similarity between data points y i and y j , other similarity measures can also be used to evaluate the similarity . For the manifold assumption, i.e., if two points y i and y j are close to each other, then their lowdimensional representations x i and x j are close as well, a natural choice is to minimize the following manifold regularization T, defined by The matrix L is usually called graph Laplacian He et al. 2011). It is apparent that minimizing the manifold regularization T imposes smoothness of the representation coefficients, or equally the prior assumption that if neighboring points y i and y j are similar (with a relatively bigger weight w ij ), their low-dimensional representations x i and x j should be very close. Therefore, minimizing (9) is an attempt to ensure the manifold assumption.
By incorporating the above manifold regularization T into the CSR, we have the MCSR problem as where MR > 0 is a manifold regularization parameter. The problem (10) has the following equivalent form where ι R + is an indicator function, defined by

Application of ADMM to MCSR
In this subsection, we propose to apply the alternating direction method of multipliers (ADMM) (Afonso et al. 2011;Yang and Zhang 2011) method to the MCSR problem (10).
The ADMM method has recently attracted more attention because it can decouple the variables, and it is usually used to solve the problems of a convex, non-smooth objective function with structured linear constraints. Consider the following structured optimization problem with linear constraints: where both f (x) and g(v) are convex functions, G is a known matrix with full column rank. For this problem, the augmented Lagrangian function is given by where α is the Lagrange multipliers, µ > 0 is a penalty parameter, and d = −α/µ. ADMM alternately minimizes L (x, v, α) with respect to x and v in a Gauss-Seidel manner. The general procedures of ADMM are summarized in Algorithm 1.
The convergence of ADMM is guaranteed by the following theorem given in Eckstein and Bertsekas (1992) (see Theorem 8).
Theorem 1 Consider problem (13) with G having full columns rank and f, g being closed, proper, convex functions. Then, for arbitrary µ > 0 and has a solution, the sequences {x t , v t , d t } generated by Algorithm 1 converges to it; otherwise, at least one of the sequences {d t } and {(x t , v t )} diverges.
According to the above framework of ADMM, and let where I is an identity matrix, and we have the corresponding augmented Lagrangian function as Then, we apply the alternating minimization idea to update the variables X, V 1 , V 2 and the Lagrange multipliers D 1 , D 2 . Given the current point X t , V t 1 , V t 2 , D t 1 , D t 2 , we get the next step of X by minimizing L with respect to X, i.e., which yields the following updating rule To update V 1 and V 2 , we have the augmented Lagrangian subproblems Before solving V 1 , we first define the well-known vector-soft thresholding (VST) operator V τ of a matrix Q = [q T 1 , q T 2 , . . . , q T m ] T ∈ R m×n as where the ith row of Q * is q * i , defined by By the VST operator, it is easy to get the updating rule of V 1 as As for V 2 , we have The Lagrangian multipliers D 1 , D 2 can be updated as Finally, the proposed manifold regularized collaborative sparse unmixing via ADMM (MCSUnADMM) algorithm for MCSR is summarized in Algorithm 2.
The convergence of Algorithm 2 is guaranteed by Theorem 1, since it can be expressed as an instance of problem (13). G is a full column rank matrix, and functions f, g are closed, proper, convex. These meet the conditions in Theorem 1, and hence the convergence of Algorithm 2 is guaranteed.
. Liu et al. SpringerPlus (2016) 5:2007 Experiments In this section, we evaluate the performance of our proposed MCSUnADD algorithm, and compare it with the CLSUnSAL algorithm (Iordache et al. 2014a) for the CSR problem (8), the SUnSAL algorithm (Iordache et al. 2011) for the SR problem (6), and the total variation regularized SUnSAL (TVSUnSAL) algorithm (Iordache et al. 2012). Experiments are carried out on two simulated hyperspectral data sets and one real hyperspectral data set, and the accuracy assessment of all the experiments is made by computing the signal to reconstruction error (SRE) (Iordache et al. 2011), defined by where E(x) is the expectation of x and x is the corresponding reconstructed abundance vector. The higher the SRE, the better the reconstruction. In our experiments, the regularization parameters SR , CSR and MR are selected from to obtain the best unmixing accuracy. As for stopping criterion, we set as the maximum iteration number and the Frobenius norm of reconstruction errors of the pixels, i.e. �Y − AX� 2 F , less than τ = 1e −4 .

Synthetic image experiments
Our experiments are first performed on two synthetic images, which were generated from the USGS library A ∈ R 224×445 that contains 445 materials with each having 224 spectral bands that uniformly distribute in the interval 0.4-2.5 μm.
Synthetic Image 1 (SI-1) This synthetic image of size 80 × 80 is produced by a hyperspectral imagery synthesis tools, which is developed by the computational intelligence group of the Basque Country University. Five endmembers as shown in Fig. 1a are randomly selected from the spectral library A to generate the synthetic image, as shown in Fig. 1b. In this image, the fractional abundances are created by the Gaussian fields method with spheric type (refer to Boris 1999 for more details) to model the natural scene such as the distribution of forest, and satisfy the LMM, together with the nonnegative and sum-to-one constraints.
Synthetic Image 2 (SI-2) This synthetic image is to simulate a scene with land covers arranged in discrete patches, and was designed in Miao and Qi (2007). In the image, seven endmembers, randomly selected from the spectral library A and as shown in Fig. 2a, are used to produce a hyperspectral image with 100 × 100 pixels. The image is first divided into several 10 pixels × 10 pixels regions with each initialized with one of the seven endmember spectra, and then a 8 pixels × 8 pixels spatial low-pass filter is used to generate the mixed pixels. In order to model the scene without pure pixels, all of the pixels whose abundances are larger than 0.8 are replaced with a mixture of all endmembers with equal abundances. By this way, the produced fractional abundances naturally satisfy the LMM with the nonnegative and sum-to-one constraints.

Results
To test the performances of algorithms influenced by the noises, zero-mean Gaussian noises are added to the above two synthetic images to achieve different signal-to-noise ratios (SNRs) of 15, 25, 35, and 45 dB. The results with SREs obtained by the MCSU-nADMM, CLSUnSAL, TVSUnSAL, and SUnSAL algorithms for these two synthetic images with different noise levels are reported in Table 1. As we can see, our proposed MCSUnADMM algorithm achieved the best unmixing accuracy than both SUnSAL and CLSUnSAL algorithms in two different simulated scenarios, and CLSUnSAL performs a little better than the SUnSAL algorithm. The TVSUnSAL algorithm has a little better performances for the SI-1 scene than that for the SI-2 scene. This may be because the two simulated scenes have different spatial characteristics (e.g., SI-1 are with a heterogeneous background while SI-2 has a relatively homogeneous background) while the TVSUnSAL algorithm aims at exploiting the spatial homogeneity to improve the unmixing accuracy. The improvement of MCSUnADMM demonstrates that the manifold regularization term is capable of enhancing the unmixing performance. Table 2 reports the times of the MCSUnADMM, CLSUnSAL, TVSUnSAL, and SUn-SAL algorithms for these two synthetic experiments with different SNRs when the stopping criterion (i.e., the maximum iteration number or the Frobenius norm of the pixel reconstruction errors, i.e. �Y − AX� 2 F , less than τ = 1e −4 ) is reached. As can been see from this table, the proposed MCSUnADMM algorithm converges faster than the other three algorithms to reach the minima of the objective function �Y − AX� 2 F because of the use of the locally geometrical structure of the hyperspectral data introduced by the manifold regularization. Note that a proper and efficient regularization introduced to the model will benefit the descending of the objective function. However, it must be pointed out that the computational complexity of solving model (10) is a little higher than solving the models (6) and (8) in theory. Table 3 gives the computational times for the SI-1 and SI-2 experiments by all of the four algorithms when setting the maximum iteration number (to be 2000) as the stopping criterion. According to this table, we can  see that the SUnSAL and CLSUnSAL algorithms need almost the same computational times, while the TVSUnSAL and the MCSUnADMM algorithms have almost the same computational times, which are higher than that of the SUnSAL and CLSUnSAL algorithms. This is consistent with the theoretical analysis of the complexity (Iordache et al. 2012(Iordache et al. , 2014a. In addition, Fig. 3 plots the evolution of the objective function (10) versus time in the SI-1 experiments with SNR = 25 dB to illustrate the convergence of the proposed MCSUnADMM algorithm. Additionally, the estimated fractional abundances obtained by the three algorithms, along with the ground-truth abundances, are shown in Figs. 4,5,6,7,8,9,10,11,12 and 13. For space consideration, only the fractional abundance maps of SI-1 with the lowest SNR of 15 dB and the fractional abundance maps of SI-2 with the highest SNR of 45 dB are reported. By visual comparisons of these fractional abundance maps, it can be seen that our proposed MCSUnADMM algorithm based on MCS model outperforms the other two algorithms for SR and CSR models and the incorporated manifold regularization can impose spatial consistency such that the spatially similar pixels have similar abundances as shown in Figs. 5,6,7,8,9 and 10. From Fig. 7, we can see that the results obtained by the TVSUnSAL algorithm exhibits more spatial homogeneity due to   the total variation regularization. However, the proposed MCSUnADMM algorithm can deal with the pixels in the transaction areas, as shown in Fig. 5, since its incorporated manifold regularization has applied different weights (9) for different pixels.

Real hyperspectral image experiments
In this section, we apply the proposed MCSUnADMM algorithm to real hyperspectral data collected by the Airborne Visible/InfRared Imaging Spectrometer (AVIRIS). The AVIRIS instrument can cover a spectral region from 0.41 to 2.45 μm in 224 bands with a 10 nm bandwidth. The ground-truth fractional abundance maps of endmember 1, 3, 5 and 7 for SI-2. a Endmember 1, b Endmember 3, c Endmember 5, d Endmember 7

Data sets
The studied real hyperspectral image is collected by AVIRIS over the Cuprite mining site, Nevada, on June 19, 1997 and is publicly available online. The Cuprite data contains some exposed minerals included in the above used USGS spectral library A, and is well understood mineralogically. The mineral map produced by USGS is shown in Fig. 14. It should be noted that we only exhibit this figure as a reference to visually compare the results obtained by different methods, since the mineral map was produced in 1995 while the Cuprite data was collected in 1997 (Iordache et al. 2014a). This scene is often used for the assessment of the abundance maps obtained by spectral unmixing algorithms (Iordache et al. 2012).

Results analysis
We conduct our experiment on a subscene with a size of 250 × 190 pixels and 224 bands. However, the bands 1-2, 104-113, 148-167, and 221-224 were removed due to water absorption and noise, and thus only a total of 188 bands were used in the experiment. In addition, the corresponding water absorption and noise bands are also removed from the spectral library A.
The fractional abundance maps of three typical minerals estimated by the MCSUn-ADMM, CLSUnSAL, and SUnSAL methods are shown in Fig. 15. By visually comparison, we can find that the abundance maps obtained by our proposed MCSUnADMM method are more sparse and have less outliers than that of the CLSUnSAL and SUnSAL methods, this may be due to the incorporated manifold regularization can consider the geometric structure of the dataset such that the performances of the MCSUnADMM method exhibit good spatial consistency.

Conclusions
This paper presents a novel sparse unmixing model, which incorporate the manifold regularization into the collaborative sparse regression. The manifold regularization is induced by a Laplacian graph, which can characterize the locally geometrical structure future work will focus on designing a strategy to adaptively set these parameters. In addition, further experiments with different scenes of real hyperspectral images are need to investigate the performances of our proposed model and method.  Fig. 15 The estimated fractional abundance maps by the MCSUnADMM (first row), CLSUnSAL (second row), and SUnSAL (third row) methods for the subscene of AVIRIS Cuprite data