Hyperspectral unmixing (HU) has become an important technique in exploiting hyperspectral data since it decomposes a mixed pixel into a collection of endmembers weighted by fractional abundances. The endmembers of a hyperspectral image (HSI) are more likely to be generated by independent sources and be mixed in a macroscopic degree before arriving at the sensor element of the imaging spectrometer as mixed spectra. Over the past few decades, many attempts have focused on imposing auxiliary regularizes on the conventional nonnegative matrix factorization (NMF) framework in order to effectively unmix these mixed spectra. As a promising step toward finding an optimum regularizer to extract endmembers, this paper presents a novel blind HU algorithm, referred to as Kurtosis-based Smooth Nonnegative Matrix Factorization (KbSNMF) which incorporates a novel regularizer based on the statistical independence of the probability density functions of endmember spectra. Imposing this regularizer on the conventional NMF framework promotes the extraction of independent endmembers while further enhancing the parts-based representation of data. Experiments conducted on diverse synthetic HSI datasets (with numerous numbers of endmembers, spectral bands, pixels, and noise levels) and three standard real HSI datasets demonstrate the validity of the proposed KbSNMF algorithm compared to several state-of-the-art NMF-based HU baselines. The proposed algorithm exhibits superior performance especially in terms of extracting endmember spectra from hyperspectral data; therefore, it could uplift the performance of recent deep learning HU methods which utilize the endmember spectra as supervisory input data for abundance extraction.

Index Terms—Hyperspectral unmixing (HU), blind source separation, kurtosis, constrained, Gaussianity, endmember indepedence, nonnegative matrix factorization (NMF).

E.M.M.B. Ekanayake is with the Department of Electrical and Computer Systems Engineering, Monash University, Clayton, Victoria 3800, Australia and also with the Office of Research and Innovation Services, Sri Lanka Technological Campus, CO 10500, Sri Lanka (email: mevan.ekanayake@monash.edu).

H.M.H.K. Weerasooriya, D.Y.L. Ranasinghe, S. Herath, G.M.R.I. Godaliyadda, H.M.V.R. Herath, and M.P.B. Ekanayake are with the Department of Electrical and Electronic Engineering, University of Peradeniya, KY 20400, Sri Lanka (email: kavingaweerasooriya@eng.pdn.ac.lk, e14273@eng.pdn.ac.lk, sanjaya.h@eng.pdn.ac.lk, roshangodd@ee.pdn.ac.lk, vijitha@eng.pdn.ac.lk; mpb.ekanayake@ee.pdn.ac.lk).

B. Rathnayake is with the Department of Electrical, Computer, and Systems Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180, USA (email: rathnb@rpi.edu).

HYPERSPECTRAL image (HSI) technology has becomea leading imaging technology in many fields including medical imaging, food quality assessment, forensic sciences, surveillance, and remote sensing [1]. However, due to the insufficient spatial resolution of spectrometers and homogeneous mixture of distinct macroscopic materials in imaging scenes, the observed reflectance spectrum at each pixel of an HSI could easily be a mixture of spectra belonging to a set of constituent members (also called endmembers). This mixing phenomenon constitutes a major concern with regard to many applications. As a remedy to this complication, various methods of hyperspectral unmixing (HU) have been implemented to extract endmember spectra along with their fractional composition (also called abundances). HU is a study of three subproblems, i.e. determining the no. of endmembers, extracting the endmember spectra, and realizing their abundances [2]. In the past, many algorithms have been introduced in order to solve the HU problem [3]–[15] and these algorithms can be categorized under three main schemes according to the basic computational approaches [16]: 1- statistical algorithms, 2- geometric algorithms and 3- sparse regression based unmxing algorithms. Statistical algorithms interpret a mixed pixel by utilizing statistical representations. These representations are commonly analytical expressions based on the probability density functions of the underlying mixed pixel spectra. Bayesian self organizing maps (BSOM) [17], independent component analysis (ICA) [18], [19], independent factor analysis (IFA) [20], dependent component analysis (DECA) [21], automated morphological endmember extraction (AMEE) [22], and spatial-spectral endmember extraction algorithm (SSEE) [23] are some of the popular statistical algorithms utilized for HU. Geometric algorithms exploit the geometric orientation of HSI data in an n-dimensional space, where n is the no. of spectral bands captured by the imaging spectrometer. Vertex component analysis (VCA) [24], minimum volume transform (MVT) [25], simplex identification via split augmented Lagrangian (SISAL) [26], optical realtime adaptive spectral identification system (ORASIS) [27], iterative error analysis (IEA) [28], and nonnegative matrix factorization (NMF) [29] are some of the geometric algorithms frequently utilized for HU. Sparse regression based approaches

utilize known libraries. The unmixing problem is formulated as a sparse linear regression problem which is based on the assumption that every feature can be linearly created by few elements extracted from known libraries [30]–[33].

In the recent past, several approaches have been introduced where deep learning (DL) is utilized for HU. More often, DLbased methods for HU do not perform blind unmixing, i.e. extract both endmember spectra and abundances. In [34], a Hopfield neural network (HNN) machine learning approach is utilized to solve the seminonnegative matrix factorization problem, which has illustrated promising performance with regard to abundance extraction when given reliable endmember spectra as supervisory input data. In [35], an artificial neural network (ANN) is utilized to inverse the pixel spectral mixture in Landsat imagery. Here, to train the network, a spectral library had been created, consisting of endmember spectra collected from the image and simulated mixed spectra. In [36], a two-staged ANN architecture has been introduced in which the first stage reduces the dimension of the input vector utilizing endmember spectra as input data. As can be seen, most of the current DL-based methods for HU utilize endmember spectra as supervisory input data in order to extract the abundances.

Originally introduced by Lee and Seung [29], NMF is a mathematical tool which is utilized to decompose a nonnegative data matrix into the product of two other nonnegative matrices of lower rank based on the optimization of a particular objective function. Since the nonnegativity criterion does not accommodate any negative elements in resultant matrices [37], which also coincides with the objective of HU. Driven by this parts-based representation of the NMF framework, NMFbased algorithms are often utilized to solve the HU problem. However, NMF is an ill-posed geometric algorithm; therefore, it does not possess a unique solution [8]. The non-convex objective functions utilized for NMF compel its solution space to be wide. Thus, many researchers have introduced novel NMF algorithms by adding different auxiliary regularizes to the conventional NMF framework in order to improve the uniqueness of its solution with respect to the HU setting. -sparsity constrained NMF (-NMF) [8], spatial group sparsity regularized NMF (SGSNMF) [38], minimum volume rank deficient NMF (Min-vol NMF) [39], manifold regularized sparse NMF [7], Double Constrained NMF [40], total variation regularized reweighted sparse NMF (TV-RSNMF) [41], subspace clustering constrained sparse NMF (SC-NMF) [42], nonsmooth NMF (nsNMF) [43], robust collaborative NMF (R-CoNMF) [44], Subspace Structure Regularized NMF (SSRNMF) [45], graph regularized NMF (GNMF) [46] and Projection-Based NMF (PNMF) [47] are some customary NMF-based baselines utilized for HU. Furthermore, A new architecture has recently emerged for blind unmixing under the premise Nonnegative Tensor factorization (NTF). In the paper, Matrix-Vector NTF for Blind Unmixing of Hyperspectral Imagery (MVNTF) [48], the authors have formalized a novel way of unmixing while preserving the spatial information by factorizing hyper spectral 3D cubes instead of unwrapped spectral datasets.

In HU, the endmembers are typically macroscopic objects

in the HSI scene, such as soil, water, vegetation, etc [2]. In a broader sense, HU attempts to find these macroscopic objects by utilizing the observations of signals that have already interacted (or mixed) with other objects in the scene before arriving at the sensing element of the imaging spectrometer. It is pragmatic to assume that the endmembers are consequences of different physical processes; hence, statistically independent1 [18]. If a particular methodology promotes maximizing the independence of endmembers, each of the endmember spectra extracted utilizing that particular method will be more independent than the mixed pixel spectra. Therefore, such a method would be a progression toward the extraction of more realistic endmember spectra belonging to independent macroscopic objects. Even though the frequently-associated abundance sum-to-one (ASC) constraint [49] in HU does not accommodate the concept of independent endmembers, algorithms such as ICA [18], [19], IFA [20], and independent innovation analysis (IIA) [50] are popular algorithms utilized in HU which consider this concept. Also, several attempts have been taken previously in order to incorporate the independence of endmembers onto the conventional NMF framework. The authors of [51] have proposed a novel initialization method based on statistical independence between NMF components. In [52], an attempt has been made to initialize NMF with a modified ICA method (modifICA-NMF). In [53], a novel effective method has been introduced unifying independent vector analysis (IVA) and NMF. Our previous work [54] discusses the suitability of utilizing the fundamental notions of kurtosis-based ICA to enhance the conventional NMF algorithm.

Inspired by the interpretable parts-based representations and simplicity of imposing auxiliary regularizes of the conventional NMF framework, and motivated by our previous work [54]–[60], this study proposes a novel regularizer to the conventional NMF framework named Average Kurtosis regularizer. Incorporating this regularizer along with an abundance smoothing mechanism, we present a novel blind HU algorithm named Kurtosis-based Smooth Nonnegative Matrix Factorization (KbSNMF) along with its two variants KbSNMF-fnorm and KbSNMF-div. The motivation of the proposed work is to promote the independence of endmembers while extracting them in accordance with the parts-based representations of the conventional NMF framework, thereby attempting to extract the most realistic endmember spectra from a given HSI. The contributions of this paper are summarized as follows:

1) Introduction of a novel regularizer for HU, based on kurtosis, which promotes the independence of endmembers of an HSI.

2) Computation of the gradient of the aforesaid regularizer w.r.t. the factors of the conventional NMF framework, and the establishment of a blind HU algorithm named KbSNMF, which effectively promotes the independence of endmembers while maintaining the smoothness of abundance maps.

We also implement and evaluate the performance of the proposed algorithm in comparison with several selected state-of-the-art NMF-based HU baselines. Experiments are conducted on diverse synthetic HSI datasets (with numerous numbers of endmembers, spectral bands, pixels, and noise levels) as well as on three standard real HSI datasets. These experiments substantiate that the proposed algorithm outperforms other state-of-the-art NMF-based blind HU algorithms in many instances, especially in extracting endmember spectra. This observation is understandable since the proposed algorithm tries to improve upon the pragmatic characteristics of the endmember spectra, rather than trying to improve upon the pragmatic characteristics of the abundance maps. Thus, in an unsupervised setting where there is the luxury of utilizing a DL-based method for abundance extraction, the proposed algorithm would provide a viable counterpart to generate endmember spectra as supervisory input data to the DL-based method.

The rest of the paper is arranged as follows. Section II provides the background related to the proposed algorithm. In Section III, the novel kurtosis-based regularizer is developed along with its derivatives. In Section IV, the novel KbSNMF algorithm is introduced. Section V discusses some key issues related to the implementation of the proposed algorithm. Section VI is devoted for experimental results and the paper is concluded in Section VII.

A. Linear Mixture Model (LMM)

LMM is the most frequently utilized model for HU and its implications had been widely discussed in many previous works [4], [6], [8], [24], [61]. This model highly depends on the assumption that the incident light waves reflect only once from the underlying macroscopic objects and are captured by the sensing element of the imaging spectrometer without being subjected to scattering. In the LMM, the spectrum at each pixel is represented as a linear combination of the endmember spectra as below,

where is the pixel spectrum, is the fractional composition occupied by the endmember in the pixel, is the spectrum of the endmember of the HSI, is an additive Gaussian noise associated with modeling errors, and r is the no. of endmembers in the HSI. All spectra are measured in reflectance values; hence, the nonnegativity in ’s and ’s. The nonnegativity constraint and the sum-to-one constraint are implied in order to guarantee that the fractional compositions representing the endmembers are nonnegative and the abundance summation equals 1 at each pixel. The LMM can be reformulated in matrix notations as below,

where is the HSI data matrix, n being the no. of spectral bands and m being the no. of pixels of the HSI, is the endmember matrix whose columns represent the spectra of each of the r endmembers, is the abundance matrix whose columns represent the fractional compositions at each of the m pixels, and is the noise matrix. This formulation casts the HU problem as a BSS problem, i.e. simultaneous extraction of the endmember spectra and their abundances at each pixel while utilizing the HSI as the input.

B. Nonnegative Matrix Factorization (NMF)

However, there are infinite no. of W, H solution pairs which satisfy the above approximation. For instance, it is possible to write for any invertible . The conventional procedure to achieve (3) is by defining an objective function which quantifies the quality of the approximation between V and WH and implementing an optimization algorithm to minimize the defined objective function w.r.t. W and H . One of the most commonly utilized objective function is the square of the Frobenius norm between V and WH as in (4).

The above expression is lower bounded by zero and distinctly vanishes if and only if V = WH. Another popular objective function is the divergence2 of V from WH as in (5).

(5) Similar to the Frobenius norm, the divergence is also lower bounded by zero and vanishes if and only if V = WH. Even though (4) and (5) functions are convex in W and H alone, they are not convex in W and H together [29]. Hence, it is not possible to analytically find global minima of these functions w.r.t. W and H. However, it is possible to find local minima utilizing numerical optimization methods. Lee and Seung [29] have proposed the below (6) and (7) multiplicative update rules to find local minima of the above (4) and (5) functions respectively.

Fig. 1: Underlying mechanism of the proposed algorithm. According to the unmixing strategy, it is discernible that every pixel is a linear combination of several independent endmembers. The proposed method promotes independence of endmembers by increasing the super-gaussianity. The algorithm concept is a derivative of the central limit theorem.

1n×mHT , H ← H ◦ WT VWHWT 1n×m (7)

Lee and Seung have further proven the convergence of both the above update rules utilizing an auxiliary function analogous to the proof of convergence of the Expectation Maximization algorithm [29].

The LMM model transforms the HU problem into the form of a conventional NMF problem. If V is the HSI data matrix X, then source matrix W is the endmember matrix A and mixing matrix H is the abundance matrix S. Thus, given X, solving the blind HU problem for A and S utilizing the conventional NMF problem can be formulated as in (8) and (9) for Frobenius norm and divergence-based objective functions respectively.

In order to solve above problems while improving the uniqueness, many previous works have incorporated additional auxiliary regularizes on A and S [8], [10], [38], [39], [44], [46].

A. Kurtosis of a Signal

Based on the value of excess kurtosis, distributions are categorized under three main types. Mesokurtic distribution is close to a Gaussian distribution; has an excess kurtosis closer to zero. Leptokurtic (also known as super-Gaussian) distribution has a higher and sharper central peak; tails are longer and flatter; has positive excess kurtosis. Platykurtic (also known as sub-Gaussian) distribution has a lower and broader central peak; tails are shorter and thinner; has negative excess kurtosis.

Central Limit Theorem (CLT) ensures that a mixture of signals is approximately Gaussian irrespective of the distributions of the underlying source signals. Even though the converse of CLT is not assured, i.e. it is not certain that any Gaussian signal is a mixture of non-Gaussian signals, in practical scenarios, Gaussian signals do consist of a mixture of non-Gaussian signals [18]. Thus, to extract the underlying source signals from a signal mixture, it is common practice in BSS to define a measure of non-Gaussianity and implement an algorithm which maximizes the defined measure as Fig. 1 illustrates. Subsequently, excess kurtosis seems to be a suitable candidate for this purpose as it is a measure of nonGaussianity. If the excess kurtosis value of a signal is close to zero, it tempts to be Gaussian and if the excess kurtosis value of a signal is away from zero, it tempts to be non-Gaussian (super- or sub-Gaussian). Since there are two types of nonGaussian distributions, it is common practice in most BSS methods to assume that source signals are of only one type. In this work, we assume the constituent spectra of an HSI to have super-Gaussian distributions. Hence, from a given HSI data matrix X, we aim to extract an endmember matrix A, whose column-wise average kurtosis is maximized, utilizing an NMF framework. Thus, we introduce a novel constrained NMF algorithm which incorporates the maximization of the average kurtosis of endmembers.

B. Average Kurtosis

Obeying the notations introduced in Section II-A, is the endmember matrix whose columns represent the spectra of each of the r endmembers of the HSI. Thus, it is possible to extract the endmember utilizing a simple matrix manipulation as below,

where is spectrum of the endmember from matrix A (or the column of matrix A) and is a column vector whose all elements are zeros except for the element which equals 1. If the kurtosis of the endmember is , it can expressed as below according to (10).

Thus, the average kurtosis through all r endmembers, K can be expressed as below utilizing (12) and (13).

Thus, it is seen that K is a function of A; therefore, it can be written as K(A). We try to maximize K so that the extracted endmembers will have a higher average kurtosis, i.e. they will be closer to super-Gaussian signals. Hence, the proposed framework would favorably influence the extraction of more realistic endmember spectra from the underlying HSI.

C. Derivative of Average Kurtosis

In order to incorporate the average kurtosis regularizer onto the conventional NMF framework, it is essential to find the gradient (or the partial derivative) of K w.r.t A and S, i.e. and . Since K is not a function of . In this section, we provide a detailed explanation on finding . Since A is the endmember matrix, we denote each of its elements by the notation , with the meaning of the reflectance value belonging to the spectral band of the endmember. Thus, the element of can be written as below implementing an element-wise derivative.

For the convenience of simplifying, we assume that each of the endmember spectra vectors have unit variance, i.e. . In order to rectify the effects of this assumption, a normalization step is carried out as discussed in Section V-B. As a result, we obtain a simplified version of as below,

where is the reflectance value belonging to the spectral band of the endmember, and is the mean reflectance of the endmember. As can be seen, is a summation of n more partial derivative terms for which the solutions can be obtained by utilizing the chain rule in calculus.

where represents a normalized version of the third central moment (skewness) of the

endmember. However, in this work we do not explore the implication of skewness within derivative of the average kurtosis. Concatenating the element-wise derivatives, we then express as the difference between two matrices as below,

Fig. 2: Effects of the smoothing parameter demonstrated on a ground truth abundance map (“Soil”) of a real HSI dataset (“Samson”): (a) = 0 (no smoothing), (b) = 0.2, (c) = 0.5, (d) = 0.8, (e) = 1 (maximum smoothing).

where and . Then, Q and P can be written as in (21) and (22) respectively for the convenience of incorporating in the NMF framework.

where and denotes a matrix whose all elements are ones. denotes the Hadamard (element-by-element) power by 3. Finally, from (20), can be written as follows.

In this section, we propose a novel blind HU algorithm which not only promotes the independence of endmembers via the kurtosis regularizer but also promotes the smoothness of the abundance maps by integrating a smoothing matrix to the conventional NMF framework. Hence, we denominate the proposed algorithm as Kurtosis-based Smooth NMF (KbSNMF). In the proceeding sections, we discuss two variants of KbSNMF depending on the objective function utilized for approximation.

A. KbSNMF-fnorm

Here we present KbSNMF based on Frobenius norm (KbSNMF-fnorm). The general optimization problem for KbSNMF-fnorm is as below.

Here, is a parameter which establishes the tradeoff between approximation error and non-Gaussianity of the endmembers rendered by K, and is a symmetric matrix called the smoothing matrix which is defined as below,

where I is the identity matrix, is a vector whose all elements are ones, and is a parameter which satisfies and controls the extent of smoothness. Enforcing smoothness onto the abundance matrix can be interpreted as Y = MS, where Y is the smoothness-enforced abundance matrix. When , hence, Y = S and no smoothing has occurred in S. As tends to become smoother and reaches the smoothest possible at . Fig. 2 demonstrates the effects of smoothing parameter on the abundance maps.

In order to find a solution for (24), we consider the objective function below.

In order to make the algorithm much simpler, the variable matrices A and S are updated in turns. In each iteration, first A is updated while S is kept constant, then, S is updated while A is kept constant. This scheme is called a block-coordinate descent approach and is widely utilized in NMFbased algorithms [66]. The updates rules can be primarily written as follows.

where denotes the Hadamard (element-by-element) product. Updating A and S directly accounts to computing the partial derivatives and , and finding suitable learning rates and .

Computing the partial derivatives of L w.r.t. A and S can be seen as two parts, i.e. partial derivatives of term and term. We refer the readers to [43] and [66] for detailed explanation of the partial derivative of . Incorporating the result in (23), we can present the partial derivatives of L as follows.

where is the scalar quantity which equals . By sub- stituting and in the original block-coordinate descent equations in (27), we can obtain the following update rules for KbSNMF-fnorm.

However, due to the subtracting terms in the gradients, the update rules (29) can enforce A and S to contain negative elements, which contradicts with the parts-based representation of the NMF framework as well as the HU setting. Thus, following a methodology similar to that proposed by Lee and Seung [29], we define data-adaptive learning rates and as below in order to ensure all positive elements in A and S at each update step.

The fraction line denotes element-by-element division. This results in the multiplicative update rules for the proposed KbSNMF-fnorm algorithm as below.

For convenience, we reconfigure the placement of matrices. Therefore, the final update rules for the proposed KbSNMFfnorm algorithm will be as follows.

It can be seen that choosing the data-adaptive learning rates in the form of (30) to avoid subtraction has enforced A and S to contain nonnegative elements throughout the block-coordinate descent approach, given initial nonnegative A and S.

B. KbSNMF-div

Analogously, we present the following optimization problem for KbSNMF based on divergence (KbSNMF-div).

Following a similar procedure as in Section IV-A, the following multiplicative update rules can be derived for KbSNMFdiv algorithm,

1n×m(MS)T + γ′N[NA]◦3

where is a matrix whose all elements are one, and the other notations are the same as defined previously.

In this section, we will discuss several points related to the implementation of the proposed algorithm.

A. Initialization

Many algorithms had been designed in the past to enhance the initialization of the conventional NMF problem. In this work, we utilize the Nonnegative Double Singular Value Decomposition (NNDSVD) algorithm [67] in order to initialize the matrices A and S. NNDSVD takes the HSI X and the no. of endmembers r as the input and generates a pair of A and S matrices. The basic NNDSVD algorithm is based on two singular value decomposition (SVD) processes, first, approximating the data matrix and the second, approximating

0 200 400 600 800 1000 iteration (a)

0 200 400 600 800 1000 iteration (e)

0 200 400 600 800 1000 iteration (c)

0 200 400 600 800 1000 iteration (b)

0 200 400 600 800 1000 iteration (f)

0 200 400 600 800 1000 iteration (d)

Fig. 3: Convergence of KbSNMF: Variation of (a) objective function, (c) average excess kurtosis of the extracted endmembers, and (e) square of Frobenius norm between X and AMS, over number iteration utilizing KbSNMFfnorm algorithm. Variation of (b) objective function, (d) average excess kurtosis of the extracted endmembers, and (f) divergence of X from AMS, over number of iteration utilizing KbSNMF-div algorithm.

Fig. 4: Standard real HSI datasets: (a) Samson dataset, (b) Urban dataset, (c) Cuprite dataset.

positive sections of the resulting partial SVD factors incorporating the properties of unit rank matrices. Extensive evidence can be found to suggest that NNDSVD promotes the rapid convergence of the NMF algorithm.

B. Normalization

To avoid the complexity of computing , the endmember spectra are considered as signals of unit variance (See Section III-C), which is not always true in HU setting. In order to rectify this premise, at the beginning of each iteration of the proposed algorithm, we normalize the endmember spectra by their individual variances (See Algorithm 1: lines 5 and 16). Thus, the resulting algorithm follows the essence of projected gradient descent methods which are often utilized in signal processing applications [18].

C. Convergence

Fig. 3 demonstrates the convergence of KbSNMF over number of iterations. Here, we have fixed the parameters and at 3 and 0.4 respectively for KbSNMF-fnorm, and at 8 and 0.4 respectively for KbSNMF-div. Selection of suitable and and their effects on the unmixing performance are extensively discussed in Section VI-C1. Observing Fig. 3(a) and 3(b), it is evident that KbSNMF convergences to a local minimum w.r.t. A and S. Also our primary objective of maximizing K has been achieved, and can be clearly seen in the Fig. 3(c) and 3(d). In the meantime, as seen in Fig. 3(e) and 3(f), Frobenius norm and divergence respectively converges to local minima w.r.t. A and S which ensures the quality of approximation between X and AMS.

D. Termination

In this work, we utilize two stopping criteria, one based on the maximum no. of iterations and the other based on the rate of change in the objective function. We choose a maximum no. of iterations, and a minimum rate of change in the objective function . The algorithm is terminated either if the present iteration t reaches or if the present rate of change in the objective function C(t) falls below . Here , where L(t) is the value of the objective function at the iteration. The selection of suitable and is discussed in Section V-E.

E. Parameter Selection

Observing Fig. 3(a) and 3(b), it is evident that both variants of KbSNMF algorithm have converged to local minima by the iteration. Thus, we fix at 1000 preserving a reasonable allowance. Also it is seen that the percentage change in the objective function around the iteration is in the order of . Thus, we fix at to ensure convergence. Determining optimum control parameters and is discussed in Section VI-C1 via experiment.

Adhering to all the implementing issues discussed above, the proposed KbSNMF algorithm can be summarized as in Algorithm 1.

A. Performance Criteria

In order to evaluate the performance of the proposed KbSNMF algorithm and assess its competitiveness with the other state-of-the-art algorithms, we utilize two performance criteria, which are commonly adopted in HU performance evaluation, i.e. Spectral Angle Distance (SAD) and Root Mean Square Error (RMSE). In most of the previous literature on HU, SAD had been utilized to compare the extracted endmember spectra with the ground truth endmember spectra while RMSE had been utilized to compare the extracted abundance maps with the ground truth abundance maps. In our work , as in (35) measures the spectral angle between the ground truth endmember spectrum and the corresponding extracted endmember spectrum , in radians; , as in (36) measures the error between the ground truth abundance map and the corresponding extracted abundance map .

Unless otherwise noted, in all experiments SAD and RMSE are average values over all extracted endmember spectra and abundance maps respectively.

B. Experimental Setting

The proposed algorithm is tested on simulated as well as real hyperspectral datasets. Also, we compare the performance of our proposed algorithm with the popular state-of-the-art NMFbased HU baselines: -NMF [8], SGSNMF [38], Min-vol NMF [39], R-CoNMF [44], SSRNMF [45] and MVNTF [48]. To ensure that the evaluations are done on common grounds, we utilize the same initializing procedure and stopping criteria as mentioned in Sections V-A and V-D respectively, for all the competing algorithms except MVNTF algorithm which is initialized with random matrices.

Simulated HSI data were generated utilizing the hyperspectral imagery synthesis toolbox (HSIST)3 in order to conduct experiments. HSIST consists of the full USGS spectral library4 which contains hundreds of endmember spectra including minerals, organic and volatile compounds, vegetation, and man-made materials. The corresponding abundance maps were generated incorporating a spherical Gaussian field [68].

To assess the performance of the proposed method in real environments, we conduct experiments on real hyperspectral

3http://www.ehu.eus/ccwintco/index.php/Hyperspectral Imagery Synth esis tools for MATLAB 4https://www.usgs.gov/labs/spec-lab

Fig. 5: Variation of unmixing performance in terms of SAD with and for (a) KbSNMF-fnorm and (b) KbSNMF-div. The minimum SAD value in each 3-D surface is marked in red.

data. The Samson dataset, the Urban dataset, and the Cuprite dataset (See Fig. 4) have been widely utilized for performance evaluation and comparison in recent HU studies [3], [42], [69]. The Samson dataset’s each pixel is recorded at 156 spectral channels covering wavelengths in the range of 401-889 nm with a spectral resolution of 3.13 nm. The Urban dataset’s each pixel is recorded at 210 spectral channels originally, however, due to dense water vapor and atmospheric effects, several bands are customarily removed prior to analysis, resulting in 162 spectral bands ranging from 400-2500 nm, with a spectral resolution of 10 nm. The Urban dataset possesses several ground truth versions, here we utilize the one with five endmembers. The Cuprite dataset is the widely used benchmark dataset for HU and its each pixel is recorded at 188 spectral channels covering wavelengths in the range of 370-2480 nm.

The ground truths for all real datasets are worked out utilizing a procedure similar to that of [4] and [70]. First, the Virtual Dimensionality (VD) algorithm [71] is utilized to determine the no. of endmembers of the HSI. Second, the pixels that contain pure endmember spectra are chosen manually in accordance with the USGS mineral spectral library. Finally, the corresponding abundances are computed utilizing

Fig. 6: Endmember spectra extracted utilizing KbSNMF: “Seawater”, “Clintonite” and “Sodiumbicarbonate” respectively.

Fig. 7: Abundance maps extracted utilizing KbSNMF: Top row- ground truth abundance maps, Middle rowextracted abundance maps by KbSNMF-fnorm. Bottom row- extracted abundance maps by KbSNMF-div.

the CVX optimization Toolbox in MATLAB. Accordinglygenerated ground truths are often utilized in HU method evaluation and comparison and are readily-available5.

C. Experiments on simulated data

1) Sensitivity to control parameters: We conduct experiments to find optimum values for and for KbSNMF-fnorm and KbSNMF-div. We increase from 0 to 25 in steps of 1, increase from 0 to 1 in steps of 0.1, and evaluate the unmixing performance at each step. It is seen that SAD reaches minimum around and in Fig. 5(a) and around and in Fig. 5(b). Thus, we fix and at 3 and 0.4 respectively for KbSNMF-fnorm and at 8 and 0.4 respectively for KbSNMF-div.

2) Unmixing performance: Under this experiment, we compare the unmixing performance of KbSNMF with the other HU algorithms. Table I shows SAD values for each of the extracted endmember spectra and Table II shows RMSE values for each of the extracted abundance maps, under the different methods. It is clearly seen that the KbSNMF under its both variants dominates the other competing algorithms in terms of SAD while signifying competitive performance in terms of RMSE. Fig. 6 and 7 respectively illustrate the endmember spectra and abundance maps extracted utilizing KbSNMF along with their ground truths.

3) Robustness to noise: In this experiment, we aim to analyze how the proposed algorithm performs in noisy environments. We add zero-mean white Gaussian noise to the original noise-free simulated dataset with a predetermined signal to noise ratio (SNR) given by (37),

where x is the pixel spectrum vector, n is the noise signal vector, and E is the expectation operator. We conduct the experiment under 11 SNR levels: 0 dB, 5 dB, 10 dB, 15 dB, 20 dB, 25 dB, 30 dB, 35 dB, 40 dB, 45 dB, 50 dB and the results are illustrated in Fig. 8 in terms of SAD and RMSE.Although SSRNMF and MVNTF show high immunity to large noise in terms of SAD values, It is discrenible that KbSNMFfnorm and KbSNMF-div report the best performance showing superior performance over all competing algorithms at noise levels in the range of 15-50 dB. They also show robustness to noise up until 30 dB. In terms of RMSE, both KbSNMF-fnorm and KbSNMF-div show robustness to noise up until 20 dB and gradually deteriorate in performance thereafter. However, both KbSNMF-fnorm and KbSNMF-div outperform SSRNMF at all noise levels in terms of RMSE. The superior performance of KbSNMF-fnorm and KbSNMF-div in terms of SAD is due to the novel auxiliary regularizes on the endmember matrix and thereby attempting to extract the most realistic endmember spectra.

4) Sensitivity to number of spectral bands: Here we vary the no. of spectral bands of the endmembers and observe the unmixing performance of the algorithms. The results are

shown in Fig. 9. KbSNMF-fnorm and KbSNMF-div outperform all the competing algorithms in terms of SAD for no. of spectral bands in the range of 300-960. However, the performance of KbSNMF-fnorm and KbSNMF-div deteriorate drastically in terms of SAD for very low no. of spectral bands, i.e. around 200 spectral bands. In terms of RMSE, KbSNMFfnorm and KbSNMF-div outperforms MVNTF at high no. of spectral bands, specifically more than 180, and outperform SGSNMF at low no. of spectral bands, i.e. below 480 spectral

Fig. 8: Variation of (a) SAD and (b) RMSE with the noise level

Fig. 9: Variation of (a) SAD and (b) RMSE with the no. of spectral bands

Fig. 10: Variation of (a) SAD and (b) RMSE with the no. of endmembers

Fig. 11: Variation of (a) SAD and (b) RMSE with the no. of pixels

Fig. 12: Endmember spectra extracted utilizing Kb- SNMF of the Samson dataset: “Soil”, “Tree” and “Water” respectively

Fig. 13: Abundance maps extracted utilizing KbSNMF of the Samson dataset: Top row- ground truth abundance maps, Middle row- extracted abundance maps by KbSNMF-fnorm. Bottom row- extracted abundance maps by KbSNMF-div.

bands.

5) Sensitivity to number of endmembers: In this experiment, we vary the no. of endmembers and investigate the performance of the algorithms. The results are illustrated in

Fig. 10. All the algorithms have the tendency to deteriorate performance in terms of SAD with the no. of endmembers. KbSNMF-fnorm and KbSNMF-div outperform R-CoNMF when the no. of endmembers are low, i.e. below 7 endmembers, and outperform SGSNMF when the no. of endmembers is high, i.e. above 5 endmembers. In terms of RMSE, KbSNMFfnorm and KbSNMF-div outperform SGSNMF when the no. of endmembers is low, i.e. below 4 endmembers.

6) Sensitivity to number of pixels: Within this experiment, we illustrate how the proposed algorithm performs under simulated HSI datasets against different no. of pixels. The no. of pixels in an HSI is a major concern since it denotes the amount of statistical information in the input to the algorithm. The amount of statistical information presented to a numerical algorithm determines the tendency of an algorithm to be trapped in a local minima [72]. Fig. 11 illustrates the results in terms of SAD and RMSE. The unmixing performance of KbSNMF-fnorm and KbSNMF-div improves in terms of SAD when the no. of pixels is increased and even outperforms all competing algorithms except MVNTF and SSRNMF when the no. of pixels is very high, i.e. 6464 and 128128 pixels. In terms of RMSE, KbSNMF-fnorm and KbSNMFdiv outperform SGSNMF when the no. of pixels is very high, i.e. 6464 and 128128 pixels.

D. Experiments on real data

We compare the unmixing performance of KbSNMF with the other competing methods in terms of SAD and RMSE for the Samson and Urban datasets. But for the Cuprite dataset, only the SAD values are tabulated.

1) Samson dataset: Table III shows SAD values for each of the extracted endmember spectra and Table IV shows RMSE values for each of the extracted abundance maps, under the different methods. In terms of average SAD, MVNTF and RCoNMF outperforms all methods. However, KbSNMF-fnorm and KbSNMF-div outperform the rest of the other methods.

TABLE I Unmixing performance comparison in terms of SAD for the simulated dataset. The best performances are in bold typeface; the second best performances are italicized; and the third best performances are underlined.

TABLE II Unmixing performance comparison in terms of RMSE for the simulated dataset. The best performances are in bold typeface; the second best performances are italicized; and the third best performances are underlined.

Fig. 14: Endmember spectra extracted utilizing Kb- SNMF of the Urban dataset: “Asphalt”, “Grass”, “Tree”, “Roof” and “Dirt” respectively.

Fig. 15: Abundance maps extracted utilizing KbSNMF of the Urban dataset: Top row- ground truth abundance maps, Middle row- extracted abundance maps by KbSNMF-fnorm. Bottom row- extracted abundance maps by KbSNMF-div.

Also, KbSNMF-div reports the third best performance in terms of SAD in extracting each endmember. In terms of RMSE, Min-vol NMF and SGSNMF outperform all methods. KbSNMF-div reports the third best average performance in terms of RMSE. The endmember spectra extracted utilizing KbSNMF-fnorm and KbSNMF-div are shown in Fig. 12, and

Fig. 16: Endmember spectra extracted utilizing Kb- SNMF of the Cuprite dataset: “Alunite”, “Andradite”, “Buddingtonite”, “Dumortierite”, “Kaolinite 1”, “Kaolinite 2”, “Muscovite”, “Montmorillonite”, “Nontronite”, “Pyrope”, “Sphene” and “Chalcedony” respectively.

it can be observed that they closely follow their ground truth spectra. Also, the abundance maps extracted by KbSNMFfnorm and KbSNMF-div are shown in Fig. 13, and it is evident that KbSNMF-div has managed to accurately extract the abundance maps.

2) Urban dataset: Table V shows SAD values for each of the extracted endmember spectra under the different methods. In terms of SAD, KbSNMF-fnorm outperforms all methods and KbSNMF-div outperforms the rest of the methods except for Min-vol NMF. Also, KbSNMF-fnorm reports the best performance and KbSNMF-div reports the second best performance in extracting the spectra of the endmembers “Tree” and “Roof”. The endmember spectra extracted utilizing KbSNMFfnorm and KbSNMF-div are shown in Fig. 14, and it can be observed that they closely follow their ground truth spectra. Table VI reports RMSE values for each of the extracted abundance maps, under the different methods. In terms of RMSE, Min-vol NMF outperforms all methods, followed by R-CoNMF and MVNTF. KbSNMF-fnorm reports the best performance in extracting the spectra of the endmembers

TABLE III Unmixing performance comparison in terms of SAD for the Samson dataset. The best performances are in bold typeface; the second best performances are italicized; and the third best performances are underlined.

TABLE IV Unmixing performance comparison in terms of RMSE for the Samson dataset; The best performances are in bold typeface; the second best performances are italicized, and the third best performances are underlined.

TABLE V Unmixing performance comparison in terms of SAD for the Urban dataset. The best performances are in bold typeface; the second best performances are italicized; and the third best performances are underlined.

TABLE VI Unmixing performance comparison in terms of RMSE for the Urban dataset; The best performances are in bold typeface; the second best performances are italicized, and the third best performances are underlined.

TABLE VII Unmixing performance comparison in terms of SAD for the Cuprite dataset; The best performances are in bold typeface; the second best performances are italicized, and the third best performances are underlined.

“Grass” and “Roof’, and KbSNMF-div reports the second best performance in extracting the spectra of the endmember “Asphalt”. The abundance maps extracted utilizing KbSNMFfnorm and KbSNMF-div are shown in Fig. 15, and it can be observed that they closely follow their ground truth abundance maps.

3) Cuprite dataset: Table VII shows SAD values for each of the extracted endmember spectra under the different methods. In terms of SAD, KbSNMF-div sits at the second place while SSRNMF stands at the top. The both KbSNMF forms report compelling results as they show best performance in extracting several endmembers, i.e “Andradite”, “Kaolinite 2”, “Muscovite”, and ”Sphene”. The endmember spectra extracted utilizing KbSNMF-fnorm and KbSNMF-div are shown in Fig. 16, and it can be observed that they closely follow their ground truth spectra.

This paper proposed a blind HU algorithm called KbSNMF, which is based on incorporating the independence of endmembers to the conventional NMF framework. This was done by introducing a novel kurtosis regularizer based on the fourth central moment of a signal which signifies the statistical independence of the underlying signal. We illustrated a comprehensive derivation of the proposed algorithm in this paper along with its performance evaluation in simulated as well as real environments (diverse simulated HSI datasets and three standard real HSI datasets). We have assessed the sensitivity of the proposed algorithm to control parameters, noise levels, number of spectral bands, number of pixels, and number of endmembers of the HSI. We have also provided performance comparisons of the proposed algorithm with the state-of-the-art NMF-based blind HU baselines. Moreover, experimental results verify that dominant performance in endmember extraction can be obtained through the novel algorithm. Hence, the proposed algorithm can be effectively utilized to extract accurate endmembers which can thereafter be passed through as supervisory input data to modern DL methods.

[1] M. J. Khan, H. S. Khan, A. Yousaf, K. Khurshid, and A. Abbas, “Modern trends in hyperspectral image analysis: A review,” IEEE Access, vol. 6, pp. 14 118–14 129, 2018.

[2] N. Keshava and J. Mustard, “Spectral unmixing,” IEEE Signal Processing Magazine, 2002.

[3] J. Qin, H. Lee, J. T. Chi, L. Drumetz, J. Chanussot, Y. Lou, and A. L. Bertozzi, “Blind hyperspectral unmixing based on graph total variation regularization,” IEEE Transactions on Geoscience and Remote Sensing, pp. 1–14, 2020.

[4] S. Jia and Y. Qian, “Constrained nonnegative matrix factorization for hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 47, no. 1, pp. 161–173, 2009.

[5] S. Yang, X. Zhang, Y. Yao, S. Cheng, and L. Jiao, “Geometric non-negative matrix factorization (gnmf) for hyperspectral unmixing,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 8, no. 6, pp. 2696–2703, 2015.

[6] F. Zhu, Y. Wang, B. Fan, S. Xiang, G. Meng, and C. Pan, “Spectral unmixing via data-guided sparsity,” IEEE Transactions on Image Processing, vol. 23, no. 12, pp. 5412–5427, 2014.

[7] X. Lu, H. Wu, Y. Yuan, P. Yan, and X. Li, “Manifold regularized sparse nmf for hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 51, no. 5, pp. 2815–2826, 2013.

[8] Y. Qian, S. Jia, J. Zhou, and A. Robles-Kelly, “Hyperspectral unmixing via sparsity-constrained nonnegative matrix factorization,” IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 11, pp. 4282–4297, Nov 2011.

[9] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret, “Nonlinear unmixing of hyperspectral images using a generalized bilinear model,” IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 11, pp. 4153–4162, 2011.

[10] C. F´evotte and N. Dobigeon, “Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization,” IEEE Transactions on Image Processing, vol. 24, no. 12, pp. 4810–4819, 2015.

[11] H. Han, G. Wang, M. Wang, J. Miao, S. Guo, L. Chen, M. Zhang, and K. Guo, “Hyperspectral unmixing via nonconvex sparse and low-rank constraint,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 13, pp. 5704–5718, 2020.

[12] N. Yokoya, T. Yairi, and A. Iwasaki, “Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion,” IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 2, pp. 528–537, 2012.

[13] M.-D. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Sparse unmixing of hyperspectral data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 6, pp. 2014–2039, 2011.

[14] Y. E. Salehani and S. Gazor, “Smooth and sparse regularization for nmf hyperspectral unmixing,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 10, no. 8, pp. 3677–3692, 2017.

[15] J. Yao, D. Meng, Q. Zhao, W. Cao, and Z. Xu, “Nonconvex-sparsity and nonlocal-smoothness-based blind hyperspectral unmixing,” IEEE Transactions on Image Processing, vol. 28, no. 6, pp. 2991–3006, 2019.

[16] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 5, no. 2, pp. 354–379, 2012.

[17] Lifan Liu, Bin Wang, Liming Zhang, and Jian Qiu Zhang, “Decomposi- tion of mixed pixels using bayesian self- organizing map (bsom) neural networks,” in 2007 IEEE International Geoscience and Remote Sensing Symposium, 2007, pp. 2014–2017.

[18] J. V. Stone, Independent Component Analysis: A Tutorial Introduction. Cambridge, MA, USA: MIT Press, 2004.

[19] J. Wang and C. . Chang, “Applications of independent component analysis in endmember extraction and abundance quantification for hyperspectral imagery,” IEEE Transactions on Geoscience and Remote Sensing, vol. 44, no. 9, pp. 2601–2616, Sep. 2006.

[20] J. M. P. Nascimento and J. M. B. Dias, “Does independent component analysis play a role in unmixing hyperspectral data?” IEEE Transactions on Geoscience and Remote Sensing, vol. 43, no. 1, pp. 175–187, Jan 2005.

[21] J. M. P. Nascimento and J. M. Bioucas-Dias, “Hyperspectral unmixing algorithm via dependent component analysis,” in 2007 IEEE International Geoscience and Remote Sensing Symposium, July 2007, pp. 4033– 4036.

[22] A. Plaza, P. Martinez, R. Perez, and J. Plaza, “Spatial/spectral endmem- ber extraction by multidimensional morphological operations,” IEEE Transactions on Geoscience and Remote Sensing, vol. 40, no. 9, pp. 2025–2041, Sep. 2002.

[23] D. Rogge, B. Rivard, J. Zhang, A. Sanchez, J. Harris, and J. Feng, “Integration of spatial–spectral information for the improved extraction of endmembers,” Remote Sensing of Environment, vol. 110, no. 3, pp. 287 – 303, 2007.

[24] J. M. P. Nascimento and J. M. B. Dias, “Vertex component analysis: a fast algorithm to unmix hyperspectral data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 43, no. 4, pp. 898–910, April 2005.

[25] M. D. Craig, “Minimum-volume transforms for remotely sensed data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 32, no. 3, pp. 542–552, May 1994.

[26] J. M. Bioucas-Dias, “A variable splitting augmented lagrangian approach

to linear spectral unmixing,” in 2009 First Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing, Aug 2009, pp. 1–4.

[27] J. H. Bowles, P. J. Palmadesso, J. A. Antoniades, M. M. Baumback, and L. J. Rickard, “Use of filter vectors in hyperspectral data analysis,” in Infrared Spaceborne Remote Sensing III, M. Strojnik and B. F. Andresen, Eds., vol. 2553, International Society for Optics and Photonics. SPIE, 1995, pp. 148 – 157.

[28] N. R. A., S. K., S. T., L. J., and H. P., “Automatic endmember extraction from hyperspectral data for mineral exploration,” in Fourth International Airborne Remote Sensing Conference and Exhibition / 21st Canadian Symposium on Remote Sensing, Ottawa, Ontario, Canada, June 1999, pp. 21–24.

[29] D. Lee and H. Seung, “Algorithms for non-negative matrix factor- ization,” in Advances in Neural Information Processing Systems 13 - Proceedings of the 2000 Conference, NIPS 2000, ser. Advances in Neural Information Processing Systems. Neural information processing systems foundation, 1 2001.

[30] M.-D. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Sparse unmixing of hyperspectral data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 6, pp. 2014–2039, 2011.

[31] M.-D. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Collaborative sparse regression for hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 52, no. 1, pp. 341–354, 2014.

[32] Z. Shi, W. Tang, Z. Duren, and Z. Jiang, “Subspace matching pursuit for sparse unmixing of hyperspectral data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 52, no. 6, pp. 3256–3274, 2014.

[33] W. Tang, Z. Shi, and Y. Wu, “Regularized simultaneous forward–backward greedy algorithm for sparse unmixing of hyperspectral data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 52, no. 9, pp. 5271–5288, 2014.

[34] J. Li, X. Li, B. Huang, and L. Zhao, “Hopfield neural network approach for supervised nonlinear spectral unmixing,” IEEE Geoscience and Remote Sensing Letters, vol. 13, no. 7, pp. 1002–1006, July 2016.

[35] Z. Mitraka, F. Del Frate, and F. Carbone, “Spectral unmixing of urban landsat imagery by means of neural networks,” in 2015 Joint Urban Remote Sensing Event (JURSE), March 2015, pp. 1–4.

[36] G. A. Licciardi and F. Del Frate, “Pixel unmixing in hyperspectral data by means of neural networks,” IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 11, pp. 4163–4172, Nov 2011.

[37] D. D. Lee and H. S. Seung, “Learning the parts of objects by non-negative matrix factorization,” Nature, vol. 401, no. 6755, pp. 788–791, 1999.

[38] X. Wang, Y. Zhong, L. Zhang, and Y. Xu, “Spatial group sparsity reg- ularized nonnegative matrix factorization for hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 55, no. 11,

pp. 6287–6304, Nov 2017.

[39] V. Leplat, A. M. S. Ang, and N. Gillis, “Minimum-volume rank-deficient nonnegative matrix factorizations,” in ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 2019, pp. 3402–3406.

[40] X. Lu, H. Wu, and Y. Yuan, “Double constrained nmf for hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 52, no. 5, pp. 2746–2758, 2014.

[41] W. He, H. Zhang, and L. Zhang, “Total variation regularized reweighted sparse nonnegative matrix factorization for hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 55, no. 7, pp. 3909–3921, 2017.

[42] X. Lu, L. Dong, and Y. Yuan, “Subspace clustering constrained sparse nmf for hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 58, no. 5, pp. 3007–3019, 2020.

[43] A. Pascual-Montano, J. M. Carazo, K. Kochi, D. Lehmann, and R. D. Pascual-Marqui, “Nonsmooth nonnegative matrix factorization (nsnmf),” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 28, no. 3, pp. 403–415, March 2006.

[44] J. Li, J. M. Bioucas-Dias, A. Plaza, and L. Liu, “Robust collaborative nonnegative matrix factorization for hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 54, no. 10, pp. 6076–6090, Oct 2016.

[45] L. Zhou, X. Zhang, J. Wang, X. Bai, L. Tong, L. Zhang, J. Zhou, and E. Hancock, “Subspace structure regularized nonnegative matrix factorization for hyperspectral unmixing,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 13, pp. 4257–4270, 2020.

[46] D. Cai, X. He, J. Han, and T. S. Huang, “Graph regularized nonnegative matrix factorization for data representation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 33, no. 8, pp. 1548– 1560, Aug 2011.

[47] Y. Yuan, Y. Feng, and X. Lu, “Projection-based nmf for hyperspectral unmixing,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 8, no. 6, pp. 2632–2643, 2015.

[48] Y. Qian, F. Xiong, S. Zeng, J. Zhou, and Y. Y. Tang, “Matrix-vector nonnegative tensor factorization for blind unmixing of hyperspectral imagery,” IEEE Transactions on Geoscience and Remote Sensing, vol. 55, no. 3, pp. 1776–1792, 2017.

[49] D. C. Heinz and Chein-I-Chang, “Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery,” IEEE Transactions on Geoscience and Remote Sensing, vol. 39, no. 3, pp. 529–545, March 2001.

[50] F. Geng, Z. Shi, Z. Jiang, and J. Yin, “Independent innovation analysis for hyperspectral imagery unmixing,” in 2008 Fourth International Conference on Natural Computation, vol. 3, Oct 2008, pp. 226–230.

[51] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and H. Saruwatari, “Determined blind source separation unifying independent vector analysis and nonnegative matrix factorization,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 24, no. 9, pp. 1622–1637, 2016.

[52] D. Benachir, S. Hosseini, Y. Deville, M. S. Karoui, and A. Hameurlain, “Modified independent component analysis for initializing non-negative matrix factorization: An approach to hyperspectral image unmixing,” in 2013 IEEE 11th International Workshop of Electronics, Control, Measurement, Signals and their application to Mechatronics, June 2013, pp. 1–6.

[53] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and H. Saruwatari, “Determined blind source separation unifying independent vector analysis and nonnegative matrix factorization,” IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 24, no. 9, pp. 1626–1641, Sep. 2016.

[54] E. M. M. B. Ekanayake, B. Rathnayake, E. M. H. E. B. Ekanayake, A. R. M. A. N. Rathnayake, H. M. V. R. Herath, G. M. R. I. Godaliyadda, and M. P. B. Ekanayake, “Enhanced hyperspectral unmixing via non-negative matrix factorization incorporating the end member independence,” in IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium, July 2019, pp. 2256–2259.

[55] B. Rathnayake, E. M. M. B. Ekanayake, K. Weerakoon, G. M. R. I. Godaliyadda, M. P. B. Ekanayake, and H. M. V. R. Herath, “Graph-based blind hyperspectral unmixing via nonnegative matrix factorization,” IEEE Transactions on Geoscience and Remote Sensing, pp. 1–19, 2020.

[56] E. M. M. B. Ekanayake, E. M. H. E. B. Ekanayake, A. R. M. A. N. Rathnayake, S. S. P. Vithana, H. M. V. R. Herath, G. M. R. I. Godaliyadda, and M. P. B. Ekanayake, “A semi-supervised algorithm to map major vegetation zones using satellite hyperspectral data,” in 2018 9th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), Sep. 2018, pp. 1–5.

[57] S. Vithana, E. Ekanayake, E. Ekanayake, A. Rathnayake, G. Jayatilaka, H. Herath, G. Godaliyadda, and M. Ekanayake, “Adaptive hierarchical clustering for hyperspectral image classification: Umbrella clustering,” Journal of Spectral Imaging, vol. 8, no. 1, p. a1, 2019.

[58] E. Ekanayake, S. Vithana, E. Ekanayake, A. Rathnayake, A. Abeysekara, T. Oorloff, H. Herath, G. Godaliyadda, M. Ekanayake, and A. Senaratne, “Mapping ilmenite deposit in pulmudai, sri lanka using a hyperspectral imaging-based surface mineral mapping method,” Journal of the National Science Foundation of Sri Lanka, vol. 47, no. 3, p. 271–284, 2019.

[59] D. Y. L. Ranasinghe, H. M. S. Lakmal, H. M. H. K. Weerasooriya, E. M. M. B. Ekanayake, G. M. R. I. Godaliyadda, H. M. V. R. Herath, and M. P. B. Ekanayake, “Hyperspectral imaging based method to identify potential limestone deposits,” in 2019 14th Conference on Industrial and Information Systems (ICIIS), 2019, pp. 135–140.

[60] E. M. M. B. Ekanayake, W. G. C. Bandara, G. W. K. Prabhath, G. M. R. I. Godaliyadda, H. M. V. R. Herath, and M. P. B. Ekanayake, “Feature extraction using minor scatter directions of data to distinguish between classes with minute differences of a hyperspectral image,” in 2019 14th Conference on Industrial and Information Systems (ICIIS), 2019, pp. 130–134.

[61] M.-D. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Total variation spatial regularization for sparse hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 11, pp. 4484–4502, 2012.

[62] L. Weixiang, Z. Nanning, and Y. Qubo, “Nonnegative matrix factor- ization and its applications in pattern recognition,” Chinese Science Bulletin, vol. 51, no. 1, pp. 7–18, 2006.

[63] A. G. K. Janecek and W. N. Gansterer, Utilizing Nonnegative Matrix Factorization for Email Classification Problems. John Wiley & Sons, Ltd, 2010, ch. 4, pp. 57–80.

[64] M. W. Berry, M. Browne, A. N. Langville, V. P. Pauca, and R. J. Plemmons, “Algorithms and applications for approximate nonnegative matrix factorization,” Computational Statistics & Data Analysis, vol. 52, no. 1, pp. 155–173, 2007.

[65] Z. Yang, G. Zhou, S. Xie, S. Ding, J. Yang, and J. Zhang, “Blind spectral unmixing based on sparse nonnegative matrix factorization,” IEEE Transactions on Image Processing, vol. 20, no. 4, pp. 1112–1125, 2011.

[66] J. J. Burred. (2014) Detailed derivation of multiplicative updaterules for nmf.

[67] C. Boutsidis and E. Gallopoulos, “Svd based initialization: A head start for nonnegative matrix factorization,” Pattern Recognition, vol. 41, no. 4, pp. 1350 – 1362, 2008.

[68] C. Shi and L. Wang, “Linear spatial spectral mixture model,” IEEE Transactions on Geoscience and Remote Sensing, vol. 54, no. 6, pp. 3599–3611, June 2016.

[69] B. Palsson, M. O. Ulfarsson, and J. R. Sveinsson, “Convolutional autoen- coder for spectral-spatial hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, pp. 1–15, 2020.

[70] S. Jia and Y. Qian, “Spectral and spatial complexity-based hyperspectral unmixing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 45, no. 12, pp. 3867–3879, Dec 2007.

[71] Chein-I Chang and Qian Du, “Estimation of number of spectrally distinct signal sources in hyperspectral imagery,” IEEE Transactions on Geoscience and Remote Sensing, vol. 42, no. 3, pp. 608–619, March 2004.

[72] W. Wang, Y. Qian, and Y. Y. Tang, “Hypergraph-regularized sparse nmf for hyperspectral unmixing,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 9, no. 2, pp. 681–694, Feb 2016.

E.M.M.B. Ekanayake (Graduate Student Member, IEEE) received the B.Sc. Engineering degree in Electrical and Electronic Engineering from the University of Peradeniya, Sri Lanka, in 2019. His undergraduate research on Hyperspectral Image Processing won the Merit Award for Manamperi Award (Engineering) in 2019 which is presented to the best undergraduate Engineering research project submitted to a Sri Lankan University. He pursued his research on Hyperspectral Image Processing as a Research Assistant at the Office of Research and Innovation Services (ORIS), Sri Lanka Technological Campus, Padukka, Sri Lanka during 2019/2020. He is currently reading for his Ph.D. at the Department of Electrical and Computer Systems Engineering, Monash University, Australia. His previous works have been published in IEEETGRS and several other IEEE-GRSS conferences including WHISPERS and IGARSS. He is a Graduate Student Member of the IEEE. His research interests include hyperspectral image processing, remote sensing, machine learning, and computer vision.

H.M.H.K. Weerasooriya obtained his degree in Electrical and Electronic Engineering with a first class honours and he currently works as an Instructor in the Department of Electronic and Electrical Engineering. Currently, he is involved in the researches on hyperspectral imaging for remote sensing and agriculture applications, and he has numerous publications in IEEE conferences. His research interests include image processing, signal processing, communication, machine learning and deep learning.

D.Y.L. Ranasinghe received his B.Sc. Engineering degree in Electrical and Electronic Engineering from the University of Peradeniya, Sri Lanka, in 2020. Immediately after, he joined the School of Engineering, Sri Lanka Technological Campus, Padukka, Sri Lanka as a Research Assistant. His research interests include hyperspectral and multispectral image analysis and processing, blind source separation, and deep learning. He has numerous publications in IEEE conferences.

S. Herath received his B.Sc. Engineering degree in Electrical and Electronic Engineering from the University of Peradeniya, Sri Lanka, in 2020. Immediately after, he joined the Department of Engineering Mathematics, University of Peradeniya as a Teaching Instructor. His research interests include computer vision, image and signal processing, pattern recognition, blind source separation, and machine learning. He has numerous publications in IEEE conferences.

B. Rathnayake received his B.Sc. Engineering degree in Electrical and Electronic Engineering from the University of Peradeniya, Sri Lanka, in 2017. Immediately after, he joined the Office of Research and Innovation Services (ORIS), Sri Lanka Technological Campus, Padukka, Sri Lanka as a Research Assistant. He is currently working as a Research Assistant at Rensselaer Polytechnic Institute, USA. His research interests include hyperspectral image analysis and processing, graph signal processing, and blind source separation. His previous works have been published in IEEE-TGRS and several other IEEE-GRSS conferences including IGARSS.

G.M.R.I. Godaliyadda (Senior Member, IEEE) obtained his B.Sc. Engineering degree in Electrical and Electronic Engineering from the University of Peradeniya, Sri Lanka, in 2005, and Ph.D. from the National University of Singapore in 2011. Currently, he is attached to the University of Peradeniya, Faculty of Engineering, Department of Electrical and Electronic Engineering as a Senior Lecturer. His current research interests include image and signal processing, pattern recognition, computer vision, machine learning, smart grid, bio-medical and remote sensing applications and algorithms. He is a Senior Member of the IEEE. He is a recipient of the Sri Lanka President’s Award for Scientific Publications for 2018 and 2019. He is the recipient of multiple grants through the National Science Foundation (NSF) for research activities. His previous works have been published in IEEE-TGRS and several other IEEE-GRSS conferences including WHISPERS and IGARSS. He also has numerous publications in many other IEEE transactions, Elsevier and IET journals and is the recipient of multiple best paper awards from international conferences for his work.

H.M.V.R. Herath (Senior Member, IEEE) received the B.Sc.Eng. degree in electrical and electronic engineering with 1st class honours from the University of Peradeniya, Peradeniya, Sri Lanka, in 1998, M.Sc. degree in electrical and computer engineering with the award of academic merit from the University of Miami, USA in 2002, and Ph.D. degree in electrical engineering from the University of Paderborn, Germany in 2009. In 2009, he joined the Department of Electrical and Electronic Engineering, University of Peradeniya, as a Senior Lecturer. His current research interests include hyperspectral imaging for remote sensing, multispectral imaging for food quality assessment, Coherent optical communications and integrated electronics. Dr. Herath was a member of one of the teams that for the first time successfully demonstrated coherent optical transmission with QPSK and polarization multiplexing. He is a member of the Institution of Engineers, Sri Lanka and The Optical Society. He is a Senior Member of the IEEE. He was the General Chair of the IEEE International Conference on Industrial and Information Systems (ICIIS) 2013 held in Kandy, Sri Lanka. His previous works have been published in IEEE-TGRS and several other IEEE-GRSS conferences including WHISPERS and IGARSS. He received the paper award in the ICTer 2017 conference held in Colombo Sri Lanka. Dr. Herath is a recipient of Sri Lanka President’s Award for scientific research in 2013.

M.P.B. Ekanayake (Senior Member, IEEE) received his B.Sc. Engineering degree in Electrical and Electronic Engineering from University of Peradeniya, Sri Lanka, in 2006, and Ph.D. from Texas Tech University in 2011. Currently, he is attached to the University of Peradeniya as a Senior Lecturer. His current research interests include applications of signal processing and system modeling in remote sensing, hyperspectral imaging, and smart grid. He is a Senior Member of the IEEE. He is a recipient of the Sri Lanka President’s Award for Scientific Publications in 2018 and 2019. He has obtained several grants through the National Science Foundation (NSF) for research projects. His previous works have been published in IEEE-TGRS and several other IEEE-GRSS conferences including WHISPERS and IGARSS. He also has multiple publications in many IEEE transactions, Elsevier and IET journals and has been awarded several best paper awards in international conferences.