Clustered Multitask Nonnegative Matrix Factorization for Spectral Unmixing of Hyperspectral Data

TL;DR

This paper introduces a novel clustered multitask nonnegative matrix factorization algorithm for spectral unmixing in hyperspectral images, utilizing fuzzy c-means clustering and diffusion LMS optimization, showing improved accuracy over existing methods.

Contribution

The paper presents a new spectral unmixing algorithm based on clustered multitask NMF with fuzzy c-means clustering and diffusion LMS, enhancing unmixing performance.

Findings

Outperforms existing methods in spectral angle distance

Achieves lower reconstruction error

Demonstrates effectiveness on synthetic and real datasets

Abstract

In this paper, the new algorithm based on clustered multitask network is proposed to solve spectral unmixing problem in hyperspectral imagery. In the proposed algorithm, the clustered network is employed. Each pixel in the hyperspectral image considered as a node in this network. The nodes in the network are clustered using the fuzzy c-means clustering method. Diffusion least mean square strategy has been used to optimize the proposed cost function. To evaluate the proposed method, experiments are conducted on synthetic and real datasets. Simulation results based on spectral angle distance, abundance angle distance and reconstruction error metrics illustrate the advantage of the proposed algorithm compared with other methods.

Tables

Table 1. Table 1: The SAD and RE performance metric and their variance (in percent) of seven algorithms on AVIRIS Cuprite dataset, using VCA-FCLS initialization.

materials	VCA-FCLS	$L_{1/2}$-NMF	GLNMF	TV-RSNMF	Dist.	S. Dist.	Proposed Al.
Sphene	0.3091	0.2143	0.1913	0.1583	0.1561	0.1673	0.1574$\pm$2.34
Nontronite	0.2622	0.2518	0.1842	0.1803	0.1944	0.1743	0.1711$\pm$1.89
KaolinSmect #1	0.2498	0.1653	0.1638	0.1731	0.2370	0.1741	0.1702$\pm$2.04
Montmorillonite	0.2609	0.2318	0.2184	0.2159	0.3571	0.2103	0.2248$\pm$2.54
Chalcedony	0.1934	0.1995	0.1649	0.1588	0.1603	0.1653	0.1437$\pm$2.61
KaolinSmect #2	0.3258	0.2542	0.2594	0.2576	0.2873	0.2608	0.2596$\pm$1.87
Alunite	0.3601	0.3458	0.2841	0.2551	0.3813	0.2369	0.2417$\pm$2.09
Buddingtonite	0.2402	0.1693	0.2068	0.2034	0.2514	0.1953	0.1643$\pm$1.91
Muscovite	0.3917	0.1584	0.1471	0.1563	0.4682	0.1537	0.1575$\pm$1.54
Andradite #1	0.2851	0.3361	0.3148	0.2392	0.2132	0.2425	0.2337$\pm$2.67
Dumortierite	0.2311	0.2453	0.2632	0.2686	0.3381	0.2639	0.2519$\pm$2.98
Andradite #2	0.4492	0.3829	0.3021	0.3136	0.3711	0.2854	0.2472$\pm$3.87
rmsSAD	0.3049	0.2562	0.2317	0.2207	0.2998	0.2153	0.2064$\pm$2.79
RE	0.0055	0.0050	0.0047	0.0038	0.0052	0.0035	0.0029$\pm$0.07

Table 2. Table 2: The SAD and RE performance metrics and their variance (in percent) of five algorithms on HYDICE Urban dataset, using VCA-FCLS initialization.

materials	VCA-FCLS	$L_{1/2}$-NMF	GLNMF	TV-RSNMF	Dist.	S. Dist.	Proposed
Roof	0.4671	0.3461	0.3486	0.3327	0.3831	0.3294	0.3289$\pm$2.56
Tree	0.2711	0.1492	0.1673	0.1572	0.2052	0.1521	0.1496$\pm$3.34
Asphalt	0.3077	0.2984	0.2096	0.2054	0.2469	0.2118	0.2122$\pm$1.96
Grass	0.2089	0.1461	0.1283	0.1249	0.1344	0.1019	0.1014$\pm$2.23
rmsSAD	0.3279	0.2512	0.2291	0.2198	0.2588	0.2161	0.2155$\pm$2.47
RE	0.0134	0.0120	0.0112	0.0108	0.0131	0.0105	0.0096$\pm$0.11

Equations

$${\mathbf{y}}_k=\mathbf{A}{\mathbf{s}}_k+{\mathbf{v}}_k$$

$$\sum _{n=1}^c{s}_k(n)=1$$

$$s_k(n)\ge 0,n=1,...,c$$

$$J^{global}({\mathbf{s}}_1,{\mathbf{s}}_2,...,{\mathbf{s}}_N)=\sum _{k=1}^N\mathrm{E}\{|{\mathbf{y}}_k-\mathbf{A}{\mathbf{s}}_k{|}^2\}$$

$${\mathbf{s}}_k^i={\mathbf{s}}_k^{i-1}+\mu \sum _{k=1}^N{\mathbf{A}}^T({\mathbf{y}}_k-\mathbf{A}{\mathbf{s}}_k)$$

$$\Delta ({\mathbf{s}}_k,{\mathbf{s}}_l)=||{\mathbf{s}}_k-{\mathbf{s}}_l|{|}^2$$

$$||\mathbf{s}_{k}||_{q}=\big{(}\sum\limits_{n=1}^{c}\mathbf{s}_{k}^{q}(n)\big{)}^{1/q}$$

$$J^{local}({\mathbf{s}}_k)=\mathrm{E}\{|{\mathbf{y}}_k-\mathbf{A}{\mathbf{s}}_k{|}^2\}+\eta \sum _{l\in {\mathcal{N}}_k\cap C(k)}{\rho }_{kl}||{\mathbf{s}}_k-{\mathbf{s}}_l|{|}^2+\lambda ||{\mathbf{s}}_k|{|}_q$$

$$\lambda =\frac{1}{\sqrt{L}}\sum _l\frac{\sqrt{N}-||{\mathbf{y}}_l|{|}_1/||{\mathbf{y}}_l|{|}_2}{\sqrt{N-1}}$$

$${\rho }_{kj}=\frac{\theta ({\mathbf{y}}_k,{\mathbf{y}}_j)}{\sum _{l\in {\mathcal{N}}_k^{-}}\theta ({\mathbf{y}}_k,{\mathbf{y}}_l)}$$

$$\theta ({\mathbf{y}}_k,{\mathbf{y}}_j)=\frac{{\mathbf{y}}_k^T{\mathbf{y}}_j}{||{\mathbf{y}}_k||||{\mathbf{y}}_j||}$$

$$\mathbf{s}_{k}^{i+1}=\mathbf{s}_{k}^{i}+\mu\mathbf{A}^{T}(\mathbf{y}_{k}-\mathbf{A}\mathbf{s}_{k})-\mu\eta\sum\limits_{l\in\mathcal{N}_{k}\cap C(k)}\rho_{kl}(\mathbf{s}_{l}^{i}-\mathbf{s}_{k}^{i})-\mu\lambda\frac{\big{(}\mathbf{s}_{k}^{i}\big{)}|\mathbf{s}_{k}^{i}|^{q-2}}{||\mathbf{s}_{k}^{i}||_{q}^{q-1}}$$

$$\underset{\mathbf{S},\mathbf{A}>0}{\min }||\mathbf{Y}-\mathbf{AS}|{|}_F^2$$

$$||\mathbf{Y}-\mathbf{AS}|{|}_F^2+\eta \sum _{k=1}^N\sum _{j\in {\mathcal{N}}_k\cap C(k)}{\rho }_{kj}||{\mathbf{s}}_j-{\mathbf{s}}_k|{|}^2+\lambda \sum _{k=1}^N||{\mathbf{s}}_k|{|}_q$$

$${\mathbf{A}}^{i+1}={\mathbf{A}}^i\ast \frac{\mathbf{Y}{\mathbf{S}}^T}{\mathbf{AS}{\mathbf{S}}^T}$$

$$\mathbf{s}_{k}^{i+1}=P^{+}\Big{(}\mathbf{s}_{k}^{i}+\mu\mathbf{A}^{T}(\mathbf{y}_{k}-\mathbf{A}\mathbf{s}_{k})-\mu\eta\sum\limits_{l\in\mathcal{N}_{k}\cap C(k)}\rho_{kl}(\mathbf{s}_{l}^{i}-\mathbf{s}_{k}^{i})-\mu\lambda\frac{\big{(}\mathbf{s}_{k}^{i}\big{)}|\mathbf{s}_{k}^{i}|^{q-2}}{||\mathbf{s}_{k}^{i}||_{q}^{q-1}}\Big{)}$$

$$||J_{new}-J_{old}||<ϵ$$

$$\text{SAD}=\cos^{-1}\Big{(}\frac{\mathbf{a}^{T}\hat{\mathbf{a}}}{||\mathbf{a}||||\hat{\mathbf{a}}||}\Big{)}$$

$$\text{AAD}=\cos^{-1}\Big{(}\frac{\mathbf{s}^{T}\hat{\mathbf{s}}}{||\mathbf{s}||||\hat{\mathbf{s}}||}\Big{)}$$

$$\text{RE}=\frac{1}{n}\sum _{i=1}^n\sqrt{||{\mathbf{y}}_i-{\hat{\mathbf{y}}}_i|{|}_2^2}$$

Full text

***Roozbeh Rajabi, \linkablerajabi@qut.ac.ir **

1 Introduction

One of the noteworthy remote sensing techniques is hyperspectral imaging. Hyperspectral images provide rich spectral information, because the sensors contain hundreds of spectral channels with higher spectral resolution than multispectral cameras. For example, the images, generating from the Airborne Visible Infra-Red Imaging Spectrometer (AVIRIS), have 224 bands. One of the problems in hyperspectral data is presence of mixed pixels [1]. In the scene, pixels containing a single material are called pure pixels and otherwise they are called mixed pixels [2]. Each pixel is composed of a set of materials called endmembers. The corresponding fraction of an endmember in that pixel is named fractional abundance [3]. Spectral unmixing (SU) methods decompose a reflectance spectrum into a set of endmember spectra [4] and their abundance fractions. The most common mixing model for hyperspectral data is linear mixing model (LMM) in which it is supposed that the recorded reflectance of a particular pixel is a linear combination of its endmembers. In contrast with nonlinear (intimate) models [5], LMM features simplicity, acceptable efficiency and low computational complexity. Because of these features, there exists many works that exploit the LMM to solve unmixing problem. Examples of such works are structured sparse method [6], minimum volume simplex analysis (MVSA) [7] that is a classic unmixing algorithm based on a minimum volume simplex, improved discrete swarm intelligence [8], conantroppy maximization using alternating direction method of multipliers (ADMM) [9] and stacked nonnegative sparse autoencoders [10] that is a special case of artificial neural network (ANN) and has the ability to extract deep robust features.

Nonnegative matrix factorization (NMF) [11, 12] is one of the practical methods of spectral unmixing, which decomposes the data into two nonnegative matrices. Recently, this basic method was developed by adding constraints, such as the minimum volume constrained NMF (MVC-NMF) method [3], graph regularized NMF (GNMF) [13], NMF with local smoothness constraint (NMF-LSC) [14], multilayer NMF (MLNMF) [15], region based structured NMF [16] and NMF based framework for hyperspectral unmixing using prior knowledge (NMFupk) [17]. Sparsity is one of the constraints for improving performance of NMF algorithm that is applied to the NMF cost function using $L_{q}$ regularizers [18], or using double reweighted sparse regression and total variation [19]. Since the number of endmembers present in each mixed pixel is small in comparison with the number of total endmembers, the problem becomes sparse [20]. $L_{1/2}$-NMF unmixing algorithm is developed by applying $L_{q}$ regularization term into NMF cost function to enforces the sparsity of endmember abundances [18].

Recently, distributed strategy gained a lot of interest in many areas [21]. Diffusion least mean square (LMS) solution has been used to solve distributed problem [22]. In a distributed optimization problem, there is a network with three types of structures: 1) a single-task network, that nodes estimate a common unknown and optimum vector, 2) a multitask network, which each node estimate its own optimum vector and 3) a clustered multitask network, which contains clusters that in each of clusters there is a common optimum vector that should be estimated [21]. In the hyperspectral image, a network of pixels can be considered to model the SU problem as a distributed one [23, 24]. Here, we have used clustered multitask network to solve SU problem [25]. Using clustered multitask network, only the neighborhood information of spectrally similar mixed pixels (those are in the same cluster) will be used in the proposed cost function. In order to generate a clustered network, we used the FCM clustering method on spectral features of hyperspectral data. Then unmixing problem has been solved as a clustered multitask network using information of nodes in neighborhood and clusters in addition to sparsity constraint.

This paper is organized as follows. In section 2, we introduce the clustered multitask network and proposed unmixing algorithm. Section 3 includes evaluation criteria, experiments on synthetic and real datasets and comparison with other methods. Section 4 gives conclusions and future work.

2 Spectral Unmixing Using Clustered Multitask Network

In this section, a new method that utilizes clustering of pixels and neighborhood information is proposed. First, we will express linear mixing model in subsection 2.1, then the distributed algorithm, cost functions and optimization procedure are formulated in 2.2. Finally, the overall algorithm to solve SU problem is presented in 2.3. The proposed method to solve spectral unmixing problem has been summarized in Figure 1.

2.1 Linear Mixing Model

To solve the SU problem, we focus on a simple and applicable model named Linear Mixing Model (LMM). In this model, observations are a linear combination of endmembers and their fractional abundances in each pixel. Mathematically, this model for pixel $k$ in the hyperspectral image is described as:

[TABLE]

where $\mathbf{y}_{k}$ is an $L\times 1$ observed data vector, $\mathbf{A}$ is the $L\times c$ signature matrix, $\mathbf{s}_{k}=[s_{k}(1),s_{k}(2),..,s_{k}(c)]^{T}$ is the $c\times 1$ fractional abundance vector and $\mathbf{v}_{k}$ is assumed as a $L\times 1$ additive noise vector of $k$-th pixel of the image, when $c$, $L$ and $N$ denote the number of endmembers, bands and pixels, respectively.

In the SU problem, fractional abundance vectors have two constraints in each pixel, abundance sum to one constraint (ASC) and abundance nonnegativity constraint (ANC) [26], which are as follows, for $c$ endmembers in a scene.

[TABLE]

[TABLE]

where $s_{k}(n)$ is the fractional abundance of the $n$-th endmember in the $k$-th pixel of the image.

2.2 Distributed Cost Functions and Optimization

As explained earlier, we can consider three types of networks containing single task, multitask and clustered multitask networks. First, $N$ nodes are considered in a clustered multitask network and a optimum vector at node $k$ is estimated. A global cost function using LMS, $J^{global}(\mathbf{s}_{1},\mathbf{s}_{2},...,\mathbf{s}_{N})$, defined as follows:

[TABLE]

where $\mathrm{E}$ is the expectation operator. Then, to minimize the cost function, the following equation is obtained, using the iterative steepest-descent solution [27]:

[TABLE]

where $\mu>0$ is a step-size parameter, and the algorithm make small jumps, using an optimum value of $\mu$. This optimum value causes stability and depends on the cost function. The algorithm will diverge with a too large value of $\mu$, and will take a long time to converge with a too small value. $i$ is iteration number.

In equation (5) the neighborhood information has not been used yet. In a distributed network, information from neighboring nodes are used to improve accuracy. In this article, we utilize the squared Euclidean distance [21]:

[TABLE]

And then, the $L_{q}$ regularizer for sparsity constraint is used [18]:

[TABLE]

Note that, the solution determined from global cost function, need to have access to information over cluster of the node, but the nodes can be considered to have availability only to information of its neighbors and the nodes of in the same cluster. Thus, for solving this problem, the following local cost function is defined, using LMS and adding the (6) and (7) constraints:

[TABLE]

where the $\mathcal{N}_{k}$ shows nodes that are in the neighborhood of node $k$, the node that exists in the cluster $C(k)$. $\eta>0$ denotes a regularization parameter [21], that controls the effect of neighborhood term, $\lambda$ is a scalar value that weights the sparsity function [18], and the nonnegative coefficients $\rho_{kl}$ are normalized spectral similarity which are obtained from correlation of data vectors [21]:

[TABLE]

where $\mathbf{y}_{l}$ denotes the $l$th band of the hyperspectral image, in this equation.

[TABLE]

where $\mathcal{N}_{k}^{-}$ include neighbors of node $k$ except itself, and $\theta$ is computed as [21]:

[TABLE]

Now, minimizing the cost function of (8), using steepest-descent algorithm, results to:

[TABLE]

Hence, this recursive equation can be used to update fractional abundance vectors in the SU problem.

2.3 Proposed Algorithm

Similar to the NMF algorithm, the least mean square error should be minimized with respect to the signatures and abundances matrices, subject to the non-negativity constraint [28]. So, the following equation is denoted, using matrix notation:

[TABLE]

where $\mathbf{A}$ and $\mathbf{S}$ are the $L\times c$ signature and $c\times N$ fractional abundances matrices, respectively, and Y denotes the $L\times N$ Hyperspectral data matrix. Then, based on described equations of the sparsity constrained distributed unmixing, the neighborhood and sparsity terms are added to (13) as follows:

[TABLE]

This cost function is minimized with respect to $\mathbf{A}$, using multiplicative update rules [28], then recursive equation of signature matrix is obtained as:

[TABLE]

And the recursive equation of fractional abundances has been obtained already in accordance with (12) as follows:

[TABLE]

where $P^{+}$ operator projects vectors onto a simplex, that adopt the ASC and ANC constraints for abundance vectors [29]. To initialize the matrices, random initialization and VCA-FCLS method can be used. Since random values may obtain local optimums, VCA-FCLS initialization method has been used in this paper. Another significant point in implementation of the algorithm is stopping criteria. This approach will be stopped until the maximum number of iteration ($T$), or the following stopping criteria is reached.

[TABLE]

where $J_{new}$ and $J_{old}$ are cost function values for two consecutive iterations and $\epsilon$ has been set to $10^{-8}$ in our experiments. Now, the proposed approach is summarized in Algorithm 1.

3 Experiments and Results

In this section, for quantitative evaluation, two common performance metrics: spectral angle distance (SAD) and abundance angle distance (AAD) [3] are used. They are defined as:

[TABLE]

[TABLE]

where $\hat{\mathbf{a}}$ is the estimation of spectral signature vectors and $\hat{\mathbf{s}}$ is the estimation of fractional abundance vectors.

Another criterion that is used for quantitative evaluation is reeconstruction error (RE) that is defined as follows [10]:

[TABLE]

where $\mathbf{y}_{i}$ and $\mathbf{\hat{y}}_{i}$ are the original and reconstructed pixels and $n$ is the total number of pixels in the hyperspectral scene.

In this article, the proposed algorithm is applied on synthetic and real datasets. First, the proposed algorithm has been applied on synthetic data. to generate this dataset, six signatures of USGS library [30] have been selected randomly, using a 7$\times$7 low pass filter and containing no pure pixels [3]. Then, the zero mean Gaussian noise with 5 different levels of SNR has been added to generated data, and performance metrics have been computed by averaging 20 Monte-Carlo runs. In Figure 2, the FCM clustering of this dataset is illustrated. Then, to choose the best number of clusters in our experiments, the SAD performance metric has been evaluated, and then according to Figure 3, the best number of clusters has been set to 6, that is equal to number of endmembers. Also, values of $\mu$ and $\eta$ has been considered equal to 0.02 and 0.1, respectively [21], and $q=2$, to gain the best results.

Then the proposed algorithm and some other algorithms: VCA-FCLS, NMF, $L_{1/2}$-NMF, distributed unmixing and sparsity constrained distributed unmixing, has been applied on the generated synthetic dataset. The comparison of performance metrics of this six different methods has been shown in Figure 4 (a) and (b), where the metrics of proposed algorithm is star-dashed line and excels other methods. Figure 5 shows the original and estimated spectral signatures for 6 endmembers, when the SNR is set to $25dB$.

Another experiment has been conducted to show the performance of the proposed algorithm for different number of pixels in the synthetic hyperspectral image. Other parameters have been set like the previous experiments. The RE metric defined in (20) is used in this experiment. Figure 7 illustrates the results of this experiment. As it can be seen in this figure, the proposed algorithm performs better when the size of hyperspectral data grows. This emphasizes that with increase of the training data size, the algorithm converges to a solution more robustly.

In the next two experiments, the proposed algorithm has been applied on real datasets: AVIRIS Cuprite and HYDICE Urban. AVIRIS Cuprite dataset is a hyperspectral data captured by the AVIRIS sensor over Cuprite, Nevada. This sensor covers wavelengths from 0.4\mu$$m to 2.5\mu$$m in 224 channels [31]. 188 bands of these 224 bands are used in the experiments and the other bands (covering bands 1, 2, 104-113, 148-167, and 221-224) have been removed which are related to water-vapor absorption or low SNR bands. Figure 6 (a) illustrates a pseudo color image of this dataset.

In HYDICE Urban dataset, There are 210 bands, that covers wavelengths from 0.4\mu$$m to 2.5\mu$$m. After removing water-vapor absorption or low SNR bands (including 1-4, 76, 87, 101-111, 136-153, and 198-210), 162 bands are used in the experiments. There are 4 distinguished materials in HYDICE Urban image: asphalt, roof, tree and grass [32]. Figure 6 (b) illustrates a pseudo color image of this real dataset.

In Figure 8, the results of FCM clustering on two real datasets is illustrated. The simulation results of spectral signatures have been shown in Figure 9 and  10. Also, $SAD$ performance metric of VCA-FCLS, $L_{1/2}$-NMF, distributed unmixing, sparsity constrained distributed unmixing and proposed method on the real datasets have been compared in Table 1 and  2, the results of proposed algorithm are available in the last column and has the best $rmsSAD$ value. In Figure 11 and  12 abundance fraction maps of the proposed method has been illustrated.

4 Conclusion

This paper proposed the clustered multitask network scheme to solve SU problem. This new algorithm considered sparsity, clustering and neighborhood information. Simulation results on synthetic and real datasets illustrated preference of the proposed approach in comparison against previously published unmixing methods in terms of SAD, AAD and RE measures. Quantitatively, the proposed method improves RE metric for about 15 percent in AVIRIS Cuprite data experiment. Despite this performance improvement, the algorithm needs FCM clustering as a preprocessing stage that will increase computational complexity and run time of the algorithm. For future works, in this paper the FCM clustering method has been used, however using more efficient clustering methods to improve the clustered network can enhance the results of the proposed method.