Fast Semi-Supervised Unmixing of Hyperspectral Image by Mutual Coherence Reduction and Recursive PCA

Abstract

Dictionary pruning step is often employed prior to the sparse unmixing process to improve the performance of library aided unmixing. This paper presents a novel recursive PCA approach for dictionary pruning of linearly mixed hyperspectral data motivated by the low-rank structure of a linearly mixed hyperspectral image. Further, we propose a mutual coherence reduction method for pre-unmixing to enhance the performance of pruning. In the pruning step we, identify the actual image endmembers utilizing the low-rank constraint. We obtain an augmented version of the data by appending each image endmember and compute PCA reconstruction error, which is a convex surrogate of matrix rank. We identify the pruned library elements according to PCA reconstruction error ratio (PRER) and PCA reconstruction error difference (PRED) and employ a recursive formulation for repeated PCA computation. Our proposed formulation identifies the exact endmember set at an affordable computational requirement. Extensive simulated and real image experiments exhibit the efficacy of the proposed algorithm in terms of its accuracy, computational complexity and noise performance.

1. Introduction

Hyperspectral imaging has attained immense popularity in remote sensing community in recent years owing to its high accuracy in classification and objbfect identification from remotely sensed images. Diverse application such as environmental studies [1], agricultural studies [2,3], mineral mapping [4], surveillance employ remotely sensed hyperspectral images. Hyperspectral images record the image intensity at several bands over the electromagnetic region [5]. The inclusion of detailed spectral information about a considerable number of spectral bands increases the discriminative ability of the imaging technology leading to higher accuracy in target detection, classification and object identification [6]. Hyperspectral imaging has been found to be very useful in identifying different objects from satellite-borne images. The object identification essentially employs the spectral unmixing method, which essentially estimates the reflectance profile of the spectrally distinct materials or endmembers. The reflectance pattern obtained from an image pixel is the resultant of reflectance profile of multiple signal sources or endmembers due to the poor spatial resolution of the imaging sensors. Spectral unmixing methods necessarily estimate the reflectance pattern of endmembers present in the image and compute its fractional abundance. Traditional unmixing methods involve three stages estimation of the number of endmembers, endmember estimation and calculation of abundance of endmembers [7].

Predominantly hyperspectral unmixing can be broadly classified into two categories unsupervised and semi-supervised unmixing [8] according to the availability of spectral library. Unsupervised unmixing methods identify the endmember and abundance matrix from the data itself. Whereas, semi-supervised approach considers the spectral library as the endmember matrix and computes the abundance matrix of the library endmembers. In recent years, semi-supervised unmixing strategy [9,10] has gained prevalence as application specific spectral libraries are available due to the rapid increase in MEMS-based optics. Dictionary pruning process identifies a smaller subset of the spectral library that can represent the image.

Traditional unsupervised unmixing methods predominantly employ convex geometric, non-negative matrix factorization or independent component analysis strategy to estimate the endmembers. The convex geometry based endmember estimation approaches include vertex component analysis [11], pixel purity index [12], convex cone analysis [13], minimum volume enclosing simplex [14,15], minimum volume simplex analysis [16], iterated constrained endmember extraction [17], simplex growing algorithm (SGA) [18]. Independent component analysis (ICA) based endmember estimation methods include [19,20,21,22]. Non-negative matrix factorization approaches for estimate the endmembers as well as abundance simultaneous incorporating regularization terms like low-rank constraints, total variation. The notable NMF based unmixing methods such as Huang et al. [23], Wang et al. [24], Tsinos et al. [25], Arngren et al. [26], Jia et al. [27], Huck et al. [28], Zhang et al. [29] employ different regularization terms to constrain the solution. However, unsupervised unmixing methods can produce satisfactory performance only when some of the image pixels contain dominant endmembers.

Many semi-supervised unmixing methods aim at computing the sparse abundance matrix assuming the spectral library as the endmember matrix. Popular sparse unmixing methods like-variable splitting and augmented Lagrangian (SUnSAL) [30] employed $l_1$ sparsity term, collaborative SUnSAL algorithm [31] combined collaborative sparse regression with the sparsity promoting term, whereas, SUnSAL-TV [32] introduced a total variation regularization term in the sparse unmixing. Among the sparse unmixing methods for abundance estimation robust sparse unmixing [33,34] method incorporates a redundant regularization term to account for endmember variability, joint local abundance method [35] performs local unmixing by exploiting structural information of image, co-evolutionary approach [36] formulates a multi-objective strategy and minimize it by evolutionary algorithm. Other works such as Feng et al. [37] proposed a spatial regularization framework which employs maximum a posteriori estimation, Themelis et al. [38] introduced a hierarchical Bayesian model based sparse unmixing method, Zhang et al. [39] transform data in framelet domain and maximize the sparsity of the obtained abundance matrix, Zhu et al. [40] proposed a correntropy maximization approach for sparse unmixing. Some recent works such as Li et al. [34], Feng et al. [41], Mei et al. [42] used spatial information alongside spectral properties of the data. Since sparse unmixing consider the whole spectral library as endmember. The prevalent sparse unmixing methods mentioned above generate an abundance matrix which has lower level of sparsity.

Some library aided unmixing methods employ a pre-processing stage, which prunes the spectral library used. Prevalent dictionary pruning based unmixing methods include orthogonal matching pursuit (OMP) [43], OMP Star [44], subspace matching pursuit (SMP) [45], compressive sampling matching pursuit (CoSaMP) [46], simultaneous orthogonal matching pursuit (SOMP) [47], MUSIC-collaborative sparse regression (MUSIC-CSR) [48], robust MUSIC-dictionary aided sparse regression (RMUSIC-DANSER) [49], sparse unmixing using spectral apriori information (SUnSPI) [50], centralized collaborative unmixing [51], deblurring and sparse unmixing [52] regularized simultaneous forward–backward greedy algorithm (RSFoBa) [53], nuclear norm approach [54]. include a pruning stage. Other works such as Li et al. [55] proposes a collaborative sparse regression approach which considers the non-linearity as an outlier and employs an inexact augmented Lagrangian method to solve the optimization problem. MUSIC-CSR algorithm [48] identifies the signal subspace and its dimension by HySIME [56] in the preliminary stage. The algorithm projects each library element on the signal subspace and identifies the signal components from the resulting projection error. Robust MUSIC algorithm (RMUSIC) [48] proposes an improved noise robust version of the inversion process, which also accounts for the variability in the reflectance profile and the discrepancy in the reflectance profile between spectral library elements and the actual image endmembers. Greedy algorithms like OMP [43], OMP star [44], SOMP [47], SMP [45], CoSaMP [46] find the best matching projections of multidimensional data onto an over-complete dictionary. However the above mentioned dictionary pruning algorithms have some inherent shortcomings, which are listed below

• The size of the pruned library for algorithms like OMP [43], SMP [45], RSFoBa [53], CoSaMP [46], SUnSPI [50] is much higher compared to the actual number of endmembers.
• Some algorithms require high computational time.
• These algorithms tend to perform poorly when the mutual coherence of library is high.

Researchers have proposed several sparse inversion approaches [31] to compute abundance of the endmembers. Among the seminal works sparse unmixing method through variable splitting and augmented Lagrangian (SUnSAL) [30] employed $l_1$ sparsity term, collaborative SUnSAL algorithm [31] added a collaborative sparse regression with the sparsity promoting term, whereas, SUnSAL-TV [32] introduced a total variation regularization term in the sparse unmixing. Among the sparse unmixing methods for abundance estimation robust sparse unmixing [33] method incorporates a redundant regularization term to account for endmember variability, joint local abundance method [35] performs local unmixing by exploiting structural information of image, co-evolutionary approach [36] formulates a multi-objective strategy and minimize it by an evolutionary optimization. Other works such as Feng et al. [37] proposed a spatial regularization framework which employs maximum a posteriori estimation, Themelis et al. [38] introduced a hierarchical Bayesian model based sparse unmixing method, Zhang et al. [39] transform data in framelet domain where abundance sparsity is maximized, Zhu et al. [40] proposed a correntropy maximization approach for sparse unmixing. Some recent works such as Li et al. [34], Feng et al. [41], Mei et al. [42] used spatial information alongside spectral properties of the data.

In this paper, we propose a novel dictionary pruning approach, where we identify the optimum image endmembers employing popular PCA based dimensionality reduction. In this work, we have employed recursive PCA formulation to minimize the computational time significantly due to the repetitive computation of eigenvalue. We also include a compressive sensing based framework to reduce the mutual coherence of spectral library. The experimental results shown in the paper demonstrate that our proposed dictionary pruning is a faster and straight-forward unmixing method, which can identify the exact endmember set.

Overall the paper is organized into the following sections Section 2 presents the signal model for linear unmixing and describes the existing algorithms, Section 3 illustrates the proposed mutual coherence reduction strategy and PCA based dictionary pruning method, Section 4 presents the results obtained on simulated as well as real images, whereas, Section 5 includes the conclusion and presents future scope of the proposed work.

2. Signal Model for Linear Unmixing

According linear mixing model the spectral reflectance profile of the i-th pixel is written as

$$\begin{array}{c}\hfill x_i=a_{i}S+w_i\end{array}$$

(1)

where, $a_i$ denotes abundance of all endmembers in the i-th pixel. $S=\left[s_1,s_2,\cdots ,s_P\right]$ is the endmember matrix which contains spectral signature of P endmembers. $w_i$ represents the noise present in the i-th pixel. The whole image $X=\left[x_1,x_2,\cdots ,x_N\right]$ consisting of N pixels is represented in matrix form as

$$\begin{array}{c}\hfill X=AS+W\end{array}$$

(2)

The abundance values satisfies abundance non-negative constraint (ANC) and abundance sum to one constraint (ASC). ANC enunciate that abundance values are non-negative, whereas, ASC indicate that abundance vectors of a pixel sum to one. These constraints are expressed as

$$\begin{array}{c}\hfill \left.\begin{array}{c}{\sum }_{i=1}^P{a}_{ij}=1\left(ASC\right)\\ 0\le a_{ij}\le 1\left(ANC\right)\end{array}\right\}\end{array}$$

(3)

2.1. Semi-Supervised Unmixing

Semi-supervised unmixing algorithms consider the whole spectral library as endmember matrix and aims to estimate abundance of the spectral library using sparse inversion. Since the spectral library employed is over-complete, the obtained abundance matrix has higher levels of sparsity, which makes abundance estimation a sparse inversion problem, which represents the data as a sparse linear mixture of the library according to

$$X=MD+W$$

(4)

where, the hyperspectral image $X=\left[x_1,x_2,\cdots ,x_N\right]$ comprises of N pixels and, the spectral library $D=\left[d_1,d_2,\cdots ,d_K\right]$ comprises of reflectance pattern of K elements; $MϵR^{N\times K}$ represents abundance matrix; and $WϵR^{N\times L}$ is the noise and residual term; Sparse unmixing algorithms obtain an abundance matrix M which leads to minimum reconstruction error while maximizing sparsity and satisfying other constraints

$$\begin{array}{c}\hfill argmi{n}_M{∥X-MD∥}_2+\lambda {∥M∥}_q\end{array}$$

(5)

where, $0\le q\le 1$ Here, the first term represents reconstruction error, whereas, the second term indicates sparsity of the obtained abundance matrix.

2.1.1. Dictionary Pruning

A hyperspectral data is represented as a mixture of pruned library as

$$\begin{array}{c}\hfill X=\widehat{M}\widehat{D}+\widehat{W}\end{array}$$

(6)

The pruned library $\widehat{D}=\left[\widehat{d_1},\widehat{d_2},\cdots ,\widehat{d_R}\right]$ contains R elements; and $\widehat{M}=\left[\widehat{m_1},\widehat{m_2},\cdots ,\widehat{m_R}\right]$ is the estimated abundance matrix; The pruned library comprises selected atoms from the spectral library ($\widehat{D}\subset D$) which can represent the image in a compact formulation. Ideally, size of the pruned library be closer to the actual number of endmembers ($R\approx P$) and $R=P$ means exact match, which us the aim of ideal dictionary pruning based semi-blind unmixing algorithms.

However, these dictionary pruning algorithms have some inherent shortcomings, which are listed below

• The size of the pruned library for algorithms like OMP [43], OMP star [44], SMP [45], RSFoBa [53], CoSaMP [46], SUnSPI [50] is much higher compared to the actual number of endmembers.
• Some algorithms require high computational time.
• These algorithms tend to perform poorly when the mutual coherence of library is high.

The mutual coherence of spectral library is the maximum cosine angle distance between any two spectral library elements. The mutual coherence of spectral library is defined as

$$\begin{array}{c}\hfill \mu \left(D\right)=arg\underset{1\le i,j\le K,i\ne j}{max}\frac{\left|d_i{}^T{d}_j\right|}{{∥d_i∥}_2{∥d_j∥}_2}\end{array}$$

(7)

The value lies in the range $\left[0,1\right]$ and higher mutual coherence indicates higher similarity between multiple atoms of the spectral library. High mutual coherence leads to the identification of endmembers with similar reflectance pattern as separately. Mutual coherence reduction leads to better dictionary pruning performance.

3. Proposed Dictionary Pruning Method

In this paper, we introduce two novel dictionary pruning algorithms PCA reconstruction error difference (PRER) and PCA reconstruction error ratio (PRER). Our proposed unmixing framework comprises of four stages noise removal, estimation of the number of endmembers, dictionary pruning and abundance computation. We include an additional mutual coherence reduction stage before unmixing for improving its performance. We utilize multi-linear regression for denoising [56], Harsanyi Ferrand Chang Virtual Dimensionality (HFC-VD) [57] for estimation of the number of endmembers along with a novel method for mutual coherence reduction. The mutual coherence reduction task have not been explored in hyperspectral sparse unmixing.

3.1. Noise Removal by Multi Linear Regression

Since efficient noise removal is pertinent to spectral unmixing we employ multilinear regression [58] framework for noise removal because of its improved performance in the hyperspectral setting [56]. This method estimates the noise present in the data by using the correlation between the consecutive spectral bands. The method models the reflectance pattern of a spectral band as a linear regressive model of other spectral bands, motivated by the high correlation between the consecutive bands.

The reflectance value of all pixels in the i-th band can be represented by

$$x_{:,i}={\beta }_i{Y}_{{\sigma }_i}+{\xi }_{:,i}$$

(8)

where, $x_{:,i}$ represent reflectance profile of the i-th band; ${\beta }_i$ is the regression coefficient; $Y_{{\sigma }_i}=[x_{:,1},x_{:,2},\cdots ,x_{:,i-1},x_{:,i+1},\cdots ,x_{:,L}]$ is the reflectance of all bands except the i-th band; and ${\xi }_{:,i}$ represents noise in the i-th band. The regression coefficient is calculated by

$$\widetilde{{\beta }_i}=X{\left(Y_{{\sigma }_i}^T{Y}_{{\sigma }_i}\right)}^{-1}{Y}_{{\sigma }_i}$$

(9)

The noise in the i-th band can be estimated as

$${\xi }_{:,i}=x_{:,i}-{\beta }_i{Y}_{{\sigma }_i}$$

(10)

The noise free image at the i-th band can be obtained by

$$\widetilde{x_{:,i}}=x_{:,i}-{\sigma }_i$$

(11)

The noise free image obtained by the process leads to improved unmixing performance.

3.2. Estimation of the Number of Endmembers

State of the art algorithms for estimation of the number of endmembers include-Harsanyi-Ferrand Chang virtual dimensionality (HFC-VD) [57], hyperspectral subspace identification by minimum error (HySIME) [56], eigen GAP index [59], eigen thresholding [60], low rank subspace estimation [61], entropy estimation of eigenvalue [62], maximal orthogonal complement algorithm (MOCA) [63], high-order statistics (HOS)-HFC [64] and hyperspectral subspace identification by SURE [65], convex geometric approach GENE-CH and GENE-AH [66], Hubness phenomenon [67] etc. In this paper, we employ HFC-VD algorithm [57] for estimation of the number of endmembers due to its accuracy in the hyperspectral setting.

3.3. Mutual Coherence Reduction

Mutual coherence of a spectral library indicates the maximum similarity between any pair of library elements. The high mutual coherence of spectral library creates complications in library aided unmixing as dictionary pruning algorithms identify consider the library elements with similar reflectance profile as distinct endmembers. Identification of duplicate endmembers reduces sparsity level of the obtained abundance matrix. The mutual coherence of a spectral library of size $K\times L$ is computed by

$$\begin{array}{c}\hfill \mu \left(\mathbf{D}\right)=arg\underset{1\le i,j\le K,i\ne j}{max}\frac{\left|d_i^T{d}_j\right|}{{∥d_i∥}_2{∥d_j∥}_2}\end{array}$$

(12)

Ideally, the performance of unmixing should remain relatively unaffected by the high mutual coherence of spectral library. Although researchers have attempted to address the problem of high mutual coherence of spectral library in sparse inversion problems and compressive sensing, its effect on hyperspectral unmixing and mutual coherence reduction task has not bee carried out in hyperspectral unmixing.

In this paper, we also introduce a compressive sensing method to reduce the mutual coherence of the spectral library used. The high mutual coherence of spectral libraries creates a challenge in the library based unmixing of hyperspectral image. Mutual coherence measure identifies the maximum degree of similarity between any pair of spectral library elements. A spectral library with high mutual coherence leads to the identification of multiple spectral library elements as single endmember.

The problem of mutual coherence reduction of dictionary or library arises in sparse inversion problems in compressive sensing setting. Compressive sensing aims at obtaining the sparsest solution to the linear system

$$\begin{array}{c}\hfill x=D\alpha \end{array}$$

(13)

in terms of $L_0$ norm. Here, $x\in {\mathcal{R}}^n$ represents the measurement data, $D\in {\mathcal{R}}^{n\times p}$ is the over-complete dictionary and $\alpha \in {\mathcal{R}}^p$ indicates the sparse coefficient vector. According to compressive sensing formulation the problem is written as

$$\begin{array}{c}\hfill \underset{\alpha }{min}{∥\alpha ∥}_0,s.t.x=D\alpha \end{array}$$

(14)

The criteria for obtaining the sparsest solution of the problem [68] are displayed below

$$\begin{array}{c}\hfill {∥\alpha ∥}_0<\frac{1}{2}\left(1+\frac{1}{\mu \left(D\right)}\right)\end{array}$$

(15)

Under this condition, $\alpha$ is the sparsest solution. The low mutual coherence of dictionary facilitates the sparsest solution whereas, high mutual coherence of dictionary creates problems in pruning. Welch et al. [69] derived a theoretical bound on mutual coherence of dictionary D of size $m\times p$. According to the bound, the minimum possible mutual coherence of the library is given by

$$\begin{array}{c}\hfill \mu \left(D\right)\ge \sqrt{\frac{p-m}{m(p-1)}}\end{array}$$

(16)

Since the dictionary employed in the process (D) is fixed, the aim of mutual coherence reduction method is to estimate an optimum projection matrix K which leads to lower values of mutual coherence ($\mu \left(M\right)$).

The mutual coherence reduction method uses a random projection matrix K as the initial transformation matrix and obtains the transformed dictionary. The transformed dictionary $M=KD$ is normalized such that the rows have unit norm.

In the mutual coherence reduction method, we minimize an alternate measure of mutual coherence termed as t-averaged mutual coherence, since, computation of mutual coherence is an NP-hard problem. Hence, we propose an alternate mutual coherence measure called t-averaged mutual coherence as this is computationally more affordable. we exploit the fact that the diagonal entries of a Gram matrix contains 1, when the library elements are normalized. The t-averaged mutual coherence [70] term is calculated according to

$$\begin{array}{c}\hfill {\mu }_t\left(M\right)=\frac{{\sum }_{1\le i,k\le j,i\ne j}{\chi }_t\left(\left|g_{ij}\right|\right)\left|g_{ij}\right|}{{\sum }_{1\le i,k\le j,i\ne j}{\chi }_t\left(\left|g_{ij}\right|\right)}\end{array}$$

(17)

We aim to minimize the mutual coherence term while satisfying the properties of Gram matrix. The mutual coherence reduction task is carried out according to the following steps in the first stage, we initialize the transformation matrix K ,normalize the rows of spectral library and compute t-averaged mutual coherence of M according to (17). In the succeeding stage, we compute the Gram matrix and shrink its elements according to

$$\begin{array}{c}\hfill g_{ij}=\left\{\begin{array}{c}\gamma g_{ij},if\left|g_{ij}\right|<\sigma \\ \gamma tsign\left(g_{ij}\right),ift>\left|g_{ij}\right|\ge \gamma t\\ g_{ij},if\gamma t\ge \left|g_{ij}\right|\end{array}\right.\end{array}$$

The shrinking or thresholding operation performed by the aforementioned process makes the matrix G a full-rank matrix. Hence, we reduce the rank of the matrix G into R by applying singular value-shrinkage. Compute the square root of G according to

$$\begin{array}{c}\hfill S^{T}S=G\end{array}$$

(18)

where $S\in {\mathcal{R}}^{L\times R}$. We minimize ${\mu }_t\left(KD\right)$ while satisfying the constraint ${∥S-KD∥}_2^2\le \xi$, which indicates that S should be a close approximate of the updated library $KD$. We employ adaptive direction method of multipliers (ADMM) [71] based optimization framework, which identifies the transformation matrix that minimizes the mutual coherence ${\mu }_t\left(PD\right)$. The optimization method uses an indirect formulation for mutual coherence reduction is as follows

$$\begin{array}{c}\hfill \underset{K}{min}{\mu }_t\left(KD\right)s.t.{∥S-KD∥}_2^2\le \xi \end{array}$$

(19)

The problem is expressed according to the Lagrangian function as

$$\begin{array}{c}\hfill \underset{K}{min}{\mu }_t\left(KD\right)+\lambda {∥S-KD∥}_2^2\end{array}$$

(20)

Here, the second term limits the power of the transformed library M. ADMM formulation employs a new slack variable and assumes that

$$\begin{array}{c}\hfill f\left(K\right)={\mu }_t\left(KD\right)\\ \hfill g\left(Z\right)={∥S-ZD∥}_2^2\end{array}$$

(21)

where, $Z=K$. ADMM framework solves the sub-problem

$$\begin{array}{c}\hfill \underset{K}{min}f\left(K\right)+g\left(Z\right),s.t.K-Z=0\end{array}$$

(22)

ADMM solution updates K, Z and U according to

$$\begin{array}{c}\hfill K^{i+1}=\underset{K}{argmin}f\left(K\right)+\frac{\rho }{2}\left({∥K-Z^i+U^i∥}_2^2\right)\\ \hfill Z^{i+1}=\underset{Z}{argmin}g\left(Z\right)+\frac{\rho }{2}\left({∥K^{i+1}-Z+U^i∥}_2^2\right)\\ \hfill U^{i+1}=U^i+K^i-Z^{i+1}\end{array}$$

(23)

The transformation matrix P obtained by the process minimizes mutual coherence of the library. The algorithmic steps are clearly described in details in Algorithm 1.

Algorithm 1: Reduction of Mutual Coherence Reduction of Spectral Library

Input: Spectral library with high mutual coherence $\mathbf{D}\in {\mathcal{R}}^{\mathbf{R}\times \mathbf{L}}$

Output: Spectral library with relatively lower mutual coherence $\varphi$

Initialization: Select a random initial projection matrix $K\in {\mathcal{R}}^{R\times R}$

1: Compute the transformed library $M=KD\in {\mathcal{R}}^{R\times L}$

2: Compute the Gram matrix $G=M^{T}M$

3: Set the threshold value t

4: Compute t-averaged mutual coherence according to Equation (17)

5: while The optimization problem Equation (20) not converged do

6:  Normalize M to unit length

7:  Shrink the elements of G according to $$g_{ij}=\{\begin{array}{c}\gamma g_{ij},if\left|g_{ij}\right|<\sigma \\ \gamma tsign\left(g_{ij}\right),ift>\left|g_{ij}\right|\ge \gamma t\\ g_{ij},if\gamma t\ge \left|g_{ij}\right|\end{array}$$

8:  Obtain the square root of the Gram matrix M according to $S^{T}S=M$

9:  Apply SVD on M and reduce the rank of M to m

10: end while

3.4. Dictionary Pruning by Recursive PCA

Any hyperspectral data lives in a substantially lower dimensional subspace, since, the data arises from a latent linear mixing process. The dimension of the subspace is close to the number of signal sources or intrinsic dimensionality of the data. Accurate identification of the intrinsic dimensionality is pivotal in dictionary pruning.

We identify the lower dimensional data subspace using Principal component analysis (PCA). Different signal processing and machine learning application have employed Principal component analysis (PCA) as a tool for dimensionality reduction. However, researchers have rarely exploited explored the possibility of employing PCA for dictionary pruning. PCA identifies a low dimensional signal subspace of dimension d from the original data space (of dimension D). These d principal components correspond to the maximum variance of the data. The first principal component represents the maximum variance, and each succeeding component corresponds to the next highest variance under the constraint that it is orthogonal to the preceding components. first-d PC’s obtained are statistically uncorrelated and orthogonal to each other. Rank d PCA minimizes the least square error such that the transformed data has low rank d

$$\begin{array}{c}\hfill \underset{\widehat{X}}{min}{∥X-\widehat{X}∥}_2^{2}s.t.rank\left(\widehat{X}\right)=d\end{array}$$

(24)

Since, PCA is a data driven transformation method, both the transformed data ($\widehat{X}$) and the reconstruction error ($E_{\left(d\right)}$) depends solely on the retained dimension (d). However, the optimum reconstruction error corresponds to the numerical rank of the data.

3.4.1. Proposed PCA Reconstruction Error Ratio Criteria (PRER)

According to Craig’s unmixing criteria [72], a hyperspectral data consisting of P endmembers lies in a $P-1$-dimensional subspace obtained by PCA transformation. Transformation of the data into $P-1$-dimension leads to optimum reconstruction error and reducing the data further do not reduce the reconstruction error significantly.

We propose a dictionary pruning idea based on PCA reconstruction error ratio. In this approach, we append each library element with the data, obtain an augmented data and transform it to $P-1$-dimension. The augmented data $Y_i=\left[X;d_i\right]$ comprise of either P endmembers or $P+1$ endmembers. We identify the number of endmembers present indirectly from PCA reconstruction error obtained from $Y_i$. Let, $E_i\left(P-1\right)$ represent the reconstruction error obtained after transforming $Y_i$ into $P-1$ dimension using PCA. Intrinsic dimensionality or numerical rank of the augmented data $Y_i$ relies on the properties of the library element added $d_i$. The numerical value of the reconstruction error obtained after transforming the augmented data also $E_i\left(P-1\right)$ depends on the properties of $d_i$. If, $d_i$ is an image endmember $E_i\left(P-1\right)$ is expected to be large, whereas, if $d_i$ is not an actual image endmember $E_i\left(P-1\right)$ is lower. We propose an index called PCA reconstruction error ratio (PRER), which is expressed as

$$\begin{array}{c}\hfill R{E}_{rat}\left(i\right)=\frac{E\left(P-1\right)}{E_i\left(P-1\right)}\end{array}$$

(25)

This index PRER has considerably lower numerical value for actual image endmembers and has a higher value for the other library elements. Hence, we consider PRER as a parameter-free indirect measure to identify the image endmembers. We present the detailed implementation of PRER based pruning in Algorithm 2.

3.4.2. Proposed PCA Reconstruction Error Difference Criteria (PRED)

Whenever a particular spectral library endmember ($d_i$) is appended with the data X, the augmented data ($Y_i=[X;d_i]$) lies in either P dimensional or $P-1$-dimensional linear subspace, depending on whether the library endmember is a part of image data or not. In the first case, when the spectral library element is also an image endmember, the intrinsic dimension of the subspace is $P-1$ otherwise, the intrinsic dimensionality is P. In the first situation the reconstruction error $E_i(P-1)$ is low, in the other scenario, $E_i(P-1)$ is much higher. The difference in reconstruction error between actual data and appended data

$$\begin{array}{c}\hfill R{E}_{dif}\left(i\right)=\left|E(P-1)-E_i(P-1)\right|\end{array}$$

(26)

gives a quantitative measure which indicates whether the spectral library endmembers are also present in the image. We present the algorithmic steps for PRED based library pruning in Algorithm 3.

Algorithm 2: PCA Reconstruction Error Ratio Criteria (PRER) for Dictionary Pruning

Input: Hyperspectral image data $\mathbf{X}\in {\mathcal{R}}^{\mathbf{N}\times \mathbf{L}}$, Spectral library $\mathbf{D}\in {\mathcal{R}}^{\mathbf{K}\times \mathbf{L}}$, Number of endmembers $\mathit{P}$

Output: Index of the pruned library $\varphi$, Pruned library $\widehat{\mathbf{D}}\in {\mathcal{R}}^{\mathbf{P}\times \mathbf{L}}$

Initialization:

1: Transform the data into $P-1$-dimension by PCA. Record PCA reconstruction error $E\left(P-1\right)$.

2: for $i<K$ and $i\leftarrow i+1$ do

3:  Append each library element with the data matrix according to ${\mathbf{Y}}_{\mathbf{i}}=\left[\mathbf{X};{\mathbf{d}}_{\mathbf{i}}\right]$

4:  Calculate the reconstruction error by transforming the appended data $Y_i$ into $P-1$-dimension by PCA. Obtained reconstruction error $E_i(P-1)$

5:  Calculate PCA reconstruction error ratio criteria ${\mathit{RE}}_{\mathit{rat}}\left(i\right)=\frac{E\left(P-1\right)}{E_i(P-1)}$

6:end for

7: Consider the $P$-elements corresponding to the minimum reconstruction error ratio ${\mathit{RE}}_{\mathit{rat}}\left(i\right)$ as endmembers. Index of pruned library $\varphi$.

8: Pruned library $\widehat{\mathbf{D}}={\mathbf{D}}_{\varphi }$

9: return Index of the pruned library elements $\varphi$, pruned library $\widehat{\mathbf{D}}$

Algorithm 3: Dictionary Pruning by PCA Reconstruction Error Difference Criteria (PRED) for Dictionary Pruning

Input: Hyperspectral image data $\mathbf{X}\in {\mathcal{R}}^{\mathbf{N}\times \mathbf{L}}$, Spectral library $\mathbf{D}\in {\mathcal{R}}^{\mathbf{K}\times \mathbf{L}}$, Number of endmembers P

Output: Index of the pruned library $\varphi$, Pruned library $\widehat{\mathbf{D}}\in {\mathcal{R}}^{P\times L}$

Initialization:

Transform the data into $P-1$-dimension by PCA. PCA reconstruction error $E\left(P-1\right)$

2: for $i<K$ and $i\leftarrow i+1$ do

Append each library element with the data matrix $Y_i=\left[X;d_i\right]$

4:  Transform the data $Y_i$ into $P-1$-dimension by PCA and record the reconstruction error $E_i(P-1)$

Calculate the difference in reconstruction error $R{E}_{dif}\left(i\right)=\left|E\left(P-1\right)-E_i(P-1)\right|$

6: end for

Consider the P-elements corresponding to the minimum reconstruction error difference $R{E}_{dif}\left(i\right)$ as image endmember index $\varphi$.

8: Obtain pruned library by $\widehat{D}=D_{\varphi }$

return Index of the pruned library elements $\varphi$ and Pruned library $\widehat{D}$

3.5. Recursive Principal Component Analysis

Our proposed library pruning methods PCA reconstruction error ratio (PRER) criteria and PCA reconstruction error ratio difference (PRED) rely on repeated computation of eigenpairs of the covariance matrix of the data. We incorporate a faster formulation to estimate the covariance matrix after augmenting a spectral library element according to rank one modification. Let, the covariance matrix of the appended data $Y_i$ be denoted by $\widehat{C_i}$. This covariance matrix after appending a row can be computed from the covariance matrix of the original data using the formula

$$\begin{array}{c}\hfill \widehat{C_i}=\frac{1}{N}C+{\left(\frac{N}{N+1}\right)}^2{d}_i^T{d}_i\end{array}$$

(27)

We perform standard eigen decomposition on this modified covariance matrix, which reduces the computational runtime.

3.6. Abundance Computation

We employ SUnSAL-TV [32] algorithm for abundance computation. Since, the hyperspectral image of any natural ground scene is smooth in the spatial domain, the abundance of the endmembers obtained by the unmixing method should also inherit the smoothness. This method exploits total variation of abundance along with sparsity and reconstruction error constraints. The overall formulation of this method is

$$\begin{array}{c}\hfill \underset{\widehat{M}}{min}{∥X-\widehat{M}\widehat{D}∥}_F^2+\lambda {∥\widehat{M}∥}_{1,1}+{\lambda }_{TV}TV\left(\widehat{M}\right)subjectto{\widehat{M}}_{ij}\ge 0\end{array}$$

(28)

Here, the first term indicates reconstruction error, whereas, the second term computes $l_1$ sparsity and the final term indicate total variation. The total variation term

$$\begin{array}{c}\hfill TV\left(\widehat{M}\right)=\sum _{i,j\in ϵ}{∥{\widehat{M}}_i-{\widehat{M}}_j∥}_1\end{array}$$

(29)

essentially represents the difference in neighbourhood pixels.

4. Results

We apply our proposed unmixing methods on a large number of synthetic and real images. We vary noise level, pixel purity level, the mutual coherence of the spectral library and number of endmembers in these synthetic image experiments.

4.1. Performance Measures

We evaluate the performance of the unmixing methods on two parameters signal to reconstruction error (SRE) and the probability of detection (Pr Det).

• Signal to Reconstruction Error (SRE)
  Signal to reconstruction error (SRE) denotes the relative power of reconstructed data with respect to the actual data.

  $$SRE=10{log}_{10}\left(\frac{{∥X∥}_2}{{∥X-\widehat{X}∥}_2}\right)$$

  (30)

  where, $\widehat{X}$ is the hyperspectral data reconstructed by the unmixing or dictionary pruning algorithm. Better unmixing leads to lower reconstruction error, which in turn increases SRE.
• Probability of Detection
  The probability of detection defines the number of spectral library endmembers accurately selected according to the formula

  $$PrDet(\widehat{\Lambda },\Lambda )=Pr\left(\widehat{\Lambda }\cap \Lambda \right)=\frac{\left|\widehat{\Lambda }\cap \Lambda \right|}{\left|A\right|}$$

  (31)

  where $\Lambda$ is the indices corresponding to the actual spectral library elements present in the image and $\widehat{\Lambda }$ is the indices corresponding to the estimated spectral library elements of the pruned library.

Value of probability of detection lies in the range $0\le PrDet\le 1$. The higher value specifies close match between actual and pruned library elements and exact match is represented by $PrDet=1$.

4.2. Description of Data

4.2.1. Synthetic Image Experiments

The synthetic images used in the experiments contain random endmembers from USGS spectral library and the abundance matrix was created according to Dirichlet distribution which satisfies ASC and ANC constraint. In the data $A1$ the number of endmembers is varied as five and ten and additive white Gaussian noise is added to the data. We alter the number of image pixels and maximum abundance of any endmember in a pixel in the synthetic data $A2$. In the data $A3$ we alter mutual coherence of the library and noise level simultaneously. In all these experiments, we perform mutual coherence reduction prior to unmixing.

4.2.2. Real Image Experiments

We used HYDICE Washington dc mall image (https://engineering.purdue.edu/biehl/MultiSpec/hyperspectral.html) and HYDICE urban image (http://lesun.weebly.com/hyperspectral-data-set.html) to validate our proposed algorithms see the Figure 1. DC Mall hyperspectral image consists of 210 spectral bands which cover the electromagnetic range 400–2400 nm. We use a 188 spectral band version, which excludes the noisy and absorption bands present in the image. The ground truth results intimate that the image endmembers are covered in USGS spectral library. HyDICE urban image (http://lesun.weebly.com/hyperspectral-data-set.html) acquired by HYDICE sensor covers the electromagnetic spectrum range 400–2500 nm and comprise of 221 spectral bands. The image has a spatial size of $200\times 200$ and consists of four endmembers as per ground-truth study [73]. However, the image contains noisy bands 1–4, 76, 87, 101–111, 136–153 and 198–210. We remove these bands before processing and use a 162 band version of the data for unmixing.

4.3. Algorithms Compared

We compare our proposed unmixing method with state of the art semi-blind unmixing algorithms like MUSIC-CSR [48], RMUSIC [49], SMP [45], RSFoBa [53] and SUnSPI [50] etc. These algorithms perform pruning of spectral library before sparse inversion.

We plot the PCA reconstruction error ratio and PCA reconstruction error difference for each spectral library elements in Figure 2, which highlights that actual library endmembers lead to lower reconstruction error ratio and reconstruction error difference. Since the PRER and PRED values for the library elements are significantly lower compared to the other library elements, it is simple to identify the endmembers from these two parameters. We display the SRE comparison of data $A1$ and $A2$ on Figure 3, which highlights that PRER and PRED obtain significantly higher SRE compared to most of the methods. Figure 3a illustrates that PRER and PRED obtain relatively higher SRE for images with high levels of mixing (lower values of maximum abundance per pixel). Figure 3b suggests that our proposed PRER and PRED outperform the prevalent methods in presence of noise. However, we do not obtain satisfactory performance on extremely high noise levels. SRE values for most of the methods predictably decrease as noise level escalates. We show the abundance images corresponding to ground truth, PRER and PRED in Figure 4, Figure 5 and Figure 6 respectively. The abundance images obtained by PRER and PRED are similar to the actual ground truth abundance image. This proves the potency of our proposed unmixing. We tabulate the probability of detection on data $A1$ and data $A2$ on

Table 1. Comparing Probability of Detection (Pr Det) for Data $A1$.

nEm	Max Abun	MUSIC-CSR	RMUSIC	SMP	RSFoBa	SUnSPI	PRER	PRED
	1	1	0.6	0.8	0.1724	0.2083	1	1
1000	0.8	0.9	0.3	0.8	0.1724	0.2	1	1
	0.6	0.9	0.4	0.6	0.1613	0.1852	1	1
	1	1	0.6	0.5	0.1563	0.2381	1	1
500	0.8	0.9	0.4	0.5	0.1926	0.2083	1	1
	0.6	0.9	0.2	0.5	0.1852	0.1926	1	0.9

and

Table 2. Comparing Probability of Deterction (Pr Det) for Data $A2$.

nEm	SNR (in dB)	MUSIC-CSR	RMUSIC	SMP	RSFoBa	SUnSPI	PRER	PRED
	No Noise	1	0.6	0.8	0.1923	0.1724	1	1
	70	1	0.3	0.8	0.1613	0.1563	1	1
	50	1	0.3	0.8	0.1515	0.1563	1	1
5	30	1	0.2	0.6	0.1471	0.1429	1	1
	20	0.8333	0.2	0.4166	0.1428	0.1351	1	0.8333
	10	0.625	0.1851	0.3846	0.1388	0.1351	0.8333	0.8333
	0	0.555	0.1786	0.3571	0.1315	0.1282	0.625	0.625
	No Noise	1	0.3	0.5	0.303	0.2381	1	1
	70	1	0.2	0.5	0.2703	0.222	1	1
	50	0.8	0.3	0.5	0.2632	0.1961	0.9	0.8
10	30	0.8	0.2	0.5	0.2326	0.1887	0.9	0.9
	20	0.7	0.2	0.3	0.1923	0.1785	0.9	0.8
	10	0.6	0.1923	0.2	0.1923	0.1724	0.8	0.8
	0	0.4	0.1887	0.1887	0.1724	0.1724	0.7	0.7

respectively. The result displayed in Table 1 illustrates that our proposed algorithms obtain probability of detection equals to almost unity. The other result displayed in Table 2 shows that PRER and PRED result into a higher probability of detection in most of the situations. However, high SNR levels-20, 10 and 0 dB make it difficult to identify the exact set of image endmembers.

Table 3. Probability of Detection for varying Mutual Coherence Created by the Synthetic Data $A3$.

Mutual Coherence	MUSIC-CSR	RMUSIC	SMP	RSFoBa	SUnSPI	PRER	PRED
1	1	0.6	0.5	0.1563	0.2381	1	1
0.8	0.8	0.5	0.6	0.1926	0.2083	1	1
0.6	0.9	0.2	0.5	0.1852	0.1926	0.9	0.9

shows that PRER and PRED obtains almost unity probability of detection even for a dictionary with high mutual coherence. We present the probability of detection result for varied mutual coherence in Table 3. This result emphasizes that PRER and PRED obtains superior unmixing performance even in the presence of spectral library with high mutual coherence level.

Our proposed PRER and PRED based unmixing relies on basic operations like covariance matrix computation and eigen decomposition, it result into low computational complexity. We have employed Cupens divide and conqer algorithm [74] for eigen decomposition. This algorithm has a computational complexity of $\mathcal{O}\left(n^{2.3}\right)$. The formulation for performing rank 1 update has computational complexity of $\mathcal{O}\left(n^3\right)$. The overall complexity of the framework is $\mathcal{O}\left(n^3\right)$. The runtime comparison of PCA was reported in [75], whereas the computational complexity of robust PCA was reported in [76]. We compare the runtime performance on an $i5$ Core 2 Duo system having 8GB RAM. The runtime plot displayed in Figure 7 illustrates that SMP [45] is the fastest, closely followed by PRER or PRED. Although SMP [45] is computationally more efficient, its moderate noise performance and lower probability of detection make it unsuitable.

5. Conclusions

This paper introduces PCA as an alternative dictionary pruning method, which accurately estimates the exact spectral library endmember set if the noise level is under certain limit and the number of endmembers present in the image is accurately estimated. We incorporate a method to reduce the mutual coherence of spectral library, which improves the unmixing performance. We also present a recursive formulation for estimation of covariance matrix after rank one modification, which significantly improves the runtime performance of the proposed method.