Efficient Weighted-Adaptive Sparse Constrained Nonnegative Tensor Factorization for Hyperspectral Unmixing

Hyperspectral unmixing aims to separate pure materials and their corresponding proportions that constitute the mixed pixels of hyperspectral imagery (HSI). Recently, the matrix-vector nonnegative tensor factorization (MV-NTF) has attracted wide attention in this field due to its natural third-tensor representation of HSI. However, the NTF-based unmixing approaches are limited to the nonunique solution and long computational time. To solve these problems, we consider the low-rank and sparse priors of each abundance map simultaneously and propose a new unmixing model adopting a weighted nuclear norm and an $L_{1/2}$ norm under the MV-NTF framework. Instead of using low-rank matrix decomposition of MV-NTF, this model imposes the low-rank property on the whole abundance map, which avoids determining the rank of the abundance map in advance. Observing that each abundance map is different, we build an adaptive update mechanism to treat each low-rank and sparse constraint differently in the model. Furthermore, a new multiplicative iterative algorithm is designed to solve the proposed model. Specially, the algorithm designed for this tensor model is simplified by using the equivalence relation with nonnegative matrix factorization. Experiments demonstrate that the proposed method is effective in improving both the unmixing effect and the solving speed.


I. INTRODUCTION
H YPERSPECTRAL imagery (HSI) is 3-D image data with rich spatial and spectral information, which has been widely studied and applied in many fields [1], [2], [3], [4]. However, due to the limitation of low spatial resolution and the complex ground distribution [5], the pixels in HSI are mostly consisted of mixtures of several different materials. Therefore, hyperspectral unmixing (HU), which decomposes mixed pixels Manuscript  into several pure materials (i.e., endmembers) and their corresponding fractions (i.e., abundances), is of great significance in improving the utility of HSI. Generally, unmixing can be categorized as linear [2] or nonlinear, including bilinear [6], [7] and multilinear [8] unmixing, relied on different interpretations of forming the mixed pixels in HSI during photon propagation. Linear mixing model (LMM) assumes that photons arriving at the sensor only interact with one material so that the pixels in HSI are linear combinations of endmembers, weighted by their corresponding abundances. The structure of LMM is simple and tractable, thus a variety of unmixing methods adopting it have emerged [9], [10], [11]. One kind of LMM-based methods is in a semisupervised fashion, which needs to get the information of endmembers in advance. Geometry-based methods have been studied to extract endmembers from HSI, such as pixel purity index [12], N-FINDR [13], and vertex component analysis (VCA) [14]. Then, the extracted endmembers are use to invert the abundances by fully constrained least squares (FCLS) [15]. These methods are premised on the presence of pure pixels in HSI. Meanwhile, the performance of extracted endmembers also affects the accuracy of subsequent abundance estimation. Sparse unmixing is another kind of semisupervised method [16], [17], which requires that the spectral library covering all endmembers in HSI is available.
Instead, there is not so much information available in real HSI scenes, thus the unsupervised methods of considering unmixing as the blind source separation came into being. Typically, nonnegative matrix factorization (NMF) [18] decomposes HSI data into the product of two nonnegative matrices representing endmember and abundance, respectively. Because of its nonnegativity and clear physical meaning, NMF and its extensions have been widely used for HU [19], [20], [21]. Unfortunately, the objective functions of NMF-based methods are usually nonconvex, which tend to fall into the local optimum. For remedying this defect, an effective way is to introduce the priors of HSI as the regular constraints of NMF. The common abundance constraints include the abundance nonnegative constraint (ANC), and the abundance sum-to-one constraint (ASC). Furthermore, sparsity, the low-rank property, and smoothness of abundances have been taken into account. Sparsity results from the fact that pixels are assumed to be mixed by a subset of endmembers so that a large number of entries in the corresponding abundance matrix are zero. Some works utilize L 0 , L 1 , and L 1/2 norms under the NMF This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ framework to enhance the sparsity of abundances [22], [23], [24], [25]. Owing to the spatial correlation among the spectral signatures of pixels, the corresponding abundance vectors are linear dependent. This fact reveals a low-rank behavior in the abundance matrix, which is often characterized by minimizing the nuclear norm [26], [27]. Besides, HSI data always vary smoothly in spatial domain, so a total variation (TV) term is exploited to promote the smoothness of abundances [28]. Moreover, there are improved NMF methods by imposing the minimum volume [29], minimum dispersion [30], or endmember dissimilarity [31] constraints upon endmembers.
These regularized NMF methods consider internal properties of HSI under matrix factorization framework. Additionally, tensor is a high-dimensional extension of matrix [32]. It is a more natural representation of HSI by directly expressing it as a third-order tensor. Therefore, more and more research works consider HU as a nonnegative tensor factorization (NTF) problem [33], [34], [35], [36]. At present, there are three typical tensor factorizations: canonical polyadic decomposition [37], tucker decomposition [38], and matrix-vector nonnegative tensor factorization (MV-NTF) [39]. Among them, MV-NTF is based on block term decomposition (BTD) [40], which regards HSI as the sum of finite third-order tensors. And each component of it is generated by the outer product of a matrix (i.e., abundance) and a vector (i.e., endmember). This model satisfies the linear spectral mixing mechanism and has attracted much attention. It is shown in [39] that the unmixing performance of MV-NTF is superior to that of the advanced NMF-based methods. However, HSI data are susceptible to noise and one limitation of MV-NTF is the lack of constraints to ensure that HSI is correctly decomposed into endmembers and abundances. Thus, Xiong et al. [41] apply the TV term to abundance maps to suppress the influence of noise. Feng et al. [42] integrate sparsity, minimum volume, and nonlinearity constraints into MV-NTF. Li et al. [43] couple NMF into MV-NTF, and impose the sparse constraint on the abundance matrix. Zheng et al. [44] add the nuclear norm and the L 1 norm to enhance the low-rank property and sparsity of abundance maps, respectively. Wang et al. [45] construct an NTF-based unmixing method with an independent constraint of endmembers. All above-mentioned methods promote the unmixing performance. However, it is still necessary to explore more detailed information to establish a more effective unmixing approach.
Recently, owing to the powerful learning and data fitting ability of deep learning strategy, various state-of-the-art learning models have been adopted for hyperspectral remote sensing [46], [47], [48]. Particularly, deep learning-based unmixing methods have developed rapidly, including using autoencoder [49], convolutional neural network [50], [51], or generative adversarial network [52], etc. In addition, they also combine additional networks to extract the characteristics of endmembers and abundances to promote unmixing performance, such as EGU-Net [53] and SeCoDe [54]. Deep learning-based unmixing methods have good applicability to the actual unmixing scene through data training. In view of data dependence and running time, however, we mainly study the traditional unmixing methods here.
Although MV-NTF considers the spatial information of HSI, the prior information of endmembers and abundances has not been fully utilized. Therefore, the focus of this study is to make full use of the spatial-spectral characteristics to obtain endmembers and abundances efficiently on the premise of protecting HSI structure. In this article, we propose a new nonnegative tensor unmixing method named efficient weighted-adaptive sparse constrained nonnegative tensor factorization (EWSP-NTF). Instead of using matrix factorization to impose a strict low-rank constraint on abundance maps in MV-NTF, we adopt a weighted nuclear norm regularization to constrain the low-rank property, which avoids estimating the rank of abundance map and preserves the details of HSI. Fig. 1 shows the low-rank and sparse characteristics of each abundance map under the MV-NTF framework. It can be seen that each abundance map has different low-rank and sparse properties. Inspired by this, we introduce the low-rank and sparse regularizations with selfadaptive weights. Differing from [43] and [44], the low-rank property and sparsity of each abundance map in this article are distinctively constrained conforming to the actual HSI scene. In lineage [45], it is effective to treat the low-rank property of each abundance map differently by minimizing the weighted nuclear norm, and we further consider the sparsity of abundance maps to improve the unmixing performance. In addition, both NMF and MV-NTF conform to the linear mixing mechanism. They are convertible when the abundance matrix of NMF is reshaped into abundance maps of MV-NTF. We take this link in the optimization process. The main contributions of this work are summarized as follows.
1) Due to the prior information that each abundance map of HSI has different low-rank and sparse properties, the weighted-adaptive low-rank and sparse regularizers are incorporated into the MV-NTF framework. Unlike the low-rank matrix decomposition used in MV-NTF, we introduce the weighted nuclear norm to constrain the low-rank property of each abundance map. And we use the L 1/2 norm to promote the sparsity of each abundance map simultaneously. Moreover, the adaptive weights are taken to treat each low-rank and sparse constraint differently in our model, which is more consistent with the real scene of HSI and more inclined to obtain the correct solution. 2) We adopt an augmented multiplicative algorithm to solve the proposed optimization problem, which introduces auxiliary variables to make the variables independent and so the subproblems are easy to solve. In order to simplify the operation of tensor decomposition, we take advantage of the equivalence relation between MV-NTF and NMF. We design the optimization algorithm under NTF framework combined with NMF, which has high computational efficiency. The rest of this article is organized as follows. Section II reviews the basic knowledge of tensor, LMM mechanism, NMF, and MV-NTF models. Section III introduces our proposed unmixing model, solving and updating processes in detail. Section IV reports the experimental results on synthetic and real-world datasets. Section V gives some conclusions and the prospects of future work. In the first row, HSI Y is decomposed into the sum of the outer product of abundance maps E r and endmembers c r . In the next two rows, the low-rank property and sparsity of each abundance map are measured by the singular values and pixel values of the matrix E r , respectively. In these graphs, the black curves depict the numerical value of each measure, and the blue curves depict the cumulative probability distribution of the values. Intuitively, the points with 95% coverage information are marked, and their corresponding dimensions are given by the pink dotted lines.

A. Notations and Definitions
Some notations and definitions of tensor algebra involved in latter work are introduced here. Tensor is a multidimensional array. In this article, scalars which are zero-order tensors are denoted by lowercase letters, e.g., x. Vectors which are firstorder tensors are denoted by boldface lowercase letters, e.g., x ∈ R I . Matrices that are second-order tensors are denoted by boldface capital letters, e.g., X ∈ R I×J . Moreover, higher order tensors that have three or more orders are identified using Euler script letters, e.g., X ∈ R I 1 ×I 2 ···×I N . Specifically, an HSI can be naturally represented as a third-order tensor Y ∈ R I 1 ×I 2 ×I 3 , which consists of I 1 × I 2 pixels and I 3 spectral bands. Three further definitions related to tensors are given as follows.
Definition 1 (Mode, Fiber, and Slice): The dimension of a tensor is called mode. An N th order tensor X ∈ R I 1 ×I 2 ×···×I N has N modes. When all modes except one are fixed in a tensor, the 1-D vector is called fiber. Furthermore, when all modes except two are fixed, the obtained 2-D matrix is called slice. Fixing two modes or one mode in a third-order tensor Y ∈ R I 1 ×I 2 ×I 3 respectively, we can get the corresponding fibers, e.g., Y i 1 ,i 2 ,: , or slices, e.g., Y :,:,i 3 .
Definition 2 (Mode-n Unfolding): The mode-n unfolding of a tensor X ∈ R I 1 ×I 2 ×···×I N refers to the matricization of the tensor, which arranges the nth mode fibers to be the columns of the resulting matrix and is denoted as X (n) . For a third-order tensor Y ∈ R I 1 ×I 2 ×I 3 , there are three matrix unfolded forms, expressed as where each of them is transformed from I 1 , I 2 , and I 3 modes, respectively.

Definition 3 (Outer Product):
The outer product of two tensors A ∈ R I 1 ×I 2 ×···×I P and B ∈ R J 1 ×J 2 ×···×J Q obtains a higher order tensor F ∈ R I 1 ×I 2 ×···×I P ×J 1 ×J 2 ×···×J Q , denoted as where each entry of F is defined as the product of corresponding entries in A and B For instance, the outer product of a matrix A ∈ R I 1 ×I 2 and a vector b ∈ R J 1 is a third-order tensor, denoted as

B. LMM, NMF, and MV-NTF
The LMM assumes that each pixel y i ∈ R K×1 with K spectral bands in HSI is a linear combination of endmembers, and its expression is . .e ri , . . ., e Ri ] T ∈ R R , and n i ∈ R K denote the endmember matrix composed of the spectral signatures of R endmembers, the abundance vector containing the proportions of all endmembers in each pixel, and the noise vector, respectively. Considering y i is one of the P pixels in the matrix Y ∈ R K×P , the mixing model (1) can be rewritten as where E = [e 1 , . . ., e i , . . ., e P ] ∈ R R×P and N ∈ R K×P denote the abundance matrix containing the abundance fractions of all pixels and the noise matrix, respectively. In the physical sense of unmixing, all above components are nonnegative. Besides, the abundance fractions of each pixel in all endmembers should be summed to one. Thus, the NMF optimization problem can be described as where · F represents the Frobenius norm, and 1 ∈ R denotes an array with all ones, in which the subscript indicates its size. Unlike NMF transforming the HSI cube into 2-D matrix for processing, MV-NTF is a decomposition model that directly builds on tensors. This model regards HSI Y ∈ R ×J×K as a third-order tensor and decomposes the tensor into the sum of R component tensors. Each component tensor is defined as the outer product of a matrix (i.e., abundance) and a vector (i.e., endmember). It is formulated as where c r ∈ R K denotes the rth endmember vector; E r ∈ R I×J denotes its corresponding abundance map, which can be made up of two low-rank matrices A r ∈ R I×L r and B r ∈ R L r ×J with L r ≤ min (I, J); N ∈ R I×J×K stands for the noise tensor. Similarly, considering the ANC and ASC in the model, the MV-NTF optimization model is where δ is a tradeoff parameter reflecting the impact of the ASC. We note that MV-NTF conforms to the LMM in a tensor form. There is an equivalent numerical relationship between (2) and (4), when taking E r as a whole. The Y in (2) is just the mode-3 unfolding Y (3) of the tensor Y when P = I × J. Besides, the abundance matrix E in (2) can be reblocked and recorded as E = [ẽ 1 , . . .,ẽ r , . . .,ẽ R ] T ∈ R R×P along row, then the vector e r ∈ R P is reshaped by the abundance map E r in (4). This link will be used to simplify the solving of the tensor model proposed later.

A. Motivation
MV-NTF is a special case of BTD, and its solution is hard to be unique unless strict conditions are met [55]. Therefore, introducing the priors of data as constraints is a direct and effective way to reduce the solution space and get the unmixing results physically meaningful. There are low-rank and sparse characteristics of abundance maps under the MV-NTF framework.
We take Jasper Ridge data as an example. As shown in Fig. 1, the low-rank property and sparsity of each abundance map are reflected by counting the number of nonzero singular values and nonzero entries, respectively. It is well known that the number of nonzero singular values is the rank of the matrix. As can be seen from the marked points in the figure, the number of singular values does not exceed or approach half when the cumulative probability of singular values reaches 95%. It shows a low-rank behavior of each abundance map. Likewise, when the cumulative probability of pixel values reaches 95%, there are only less than half of nonzero entries. It shows that most entries in the abundance map are zero conversely, which indicates a sparsity of the abundance map. Furthermore, Fig. 1 shows that each abundance map corresponds to a different value when the cumulative probability distribution reaches 95%. It reflects that the low-rank property and sparsity of each abundance map are not completely the same in degree. For making full use of this detailed abundance information, we introduce the weighted-adaptive low-rank and sparse regularizers in our model. Besides, the forceful low-rank constraint of MV-NTF is replaced to protect the details of HSI.

B. Proposed EWSP-NTF Model
Based on above facts, we adopt a weighted nuclear norm and an L 1/2 norm with adaptive weights to constrain the characteristics. The proposed model is established as follows: where 1 I×J denotes an I × J matrix with all entries as one; δ, λ, and ρ are parameters of regularizations; α r and β r are self-adaptive parameters, which vary with the low-rank property and sparsity of each abundance map. The definitions of weighted nuclear norm and L 1/2 norm are given as follows.
For a matrix X ∈ R I 1 ×I 2 , the rank constraint is usually replaced by the nuclear norm that represents the sum of singular values, i.e., X * = i σ i (X), where σ i (X) is the ith singular value of X. Since the larger singular value indicates the more important information in image data, the weighted nuclear norm is introduced to treat singular values differently [56]. It is defined as where w i is the nonnegative weight coefficient of σ i (X). The L 1/2 norm is used to describe the sparsity of data [24], and it is defined as where x i 1 i 2 denotes the entry in row i 1 and column i 2 of matrix X.

C. Updating Rules
The augmented multiplicative algorithm is used to solve (6). We introduce auxiliary variables U r and V r for E r , then the equivalent form can be obtained as follows: Subsequently, constructing the Lagrange function, (9) can be reformulated as where are consistent with the expression in Section II-B; Ψ and Θ are Lagrange multipliers indicating the nonnegative variable constraints in (9), and μ > 0 is a penalty parameter. Now, the optimization problem can be separated into independent subproblems. The variables E, C, U r , and V r are updated alternately in the following.
Update E: With fixed C, U r , and V r , the suboptimization problem for E is Here, we take advantage of the equivalence relation between MV-NTF and NMF pointed out in Section II-B to unify the variables. Considering this convex subproblem (11) with Karush-Kuhn-Tucker conditions, then we have When Ψ. * E = 0, the problem (11) gets its optimal value. Hence, the update rule for E is derived as Update C: With fixed E, U r , and V r , the suboptimization problem for C is where E and C are the same as (11). Similarly, we have the linear equation where Θ. * C = 0, and the update rule for C is obtained Update U r : With fixed E, C, and V r , the suboptimization problem for U r is We solve it by using soft-thresholding operator of the singular values. Define the shrinkage operator as propose to update U r as where w r is an adaptive lowrank factor of the abundance map, and it is updated by where also set = 10 −16 . Update V r : With fixed E, C, and U r , the suboptimization problem for V r is Calculating the partial derivative of V r , we have Then the update rule of V r can be obtained as where β (t) r is an adaptive sparsity factor of the abundance map, and it is updated by with = 10 −16 . We summarize above updating rules of the proposed EWSP-NTF unmixing algorithm in Algorithm 1.

D. Computational Complexity and Implementation Issues
First, we briefly analyze the computational complexity of Algorithm 1. Instead of using the nested matrix operations in most MV-NTF-based methods, the computational cost of EWSP-NTF method is mainly contributed by matrix multiplications which are smaller in scale. As shown in Algorithm 1, each iteration includes four updating steps: for (13), the number of floating-point operations required is (2K + 11)IJR + (2K + 4IJ + 1)R 2 to update E, and for (16), the number of floating-point operations required is (2IJ + 2)KR + (2K + 2IJ)R 2 to update C; since U and V are divided into R submatrices U r and V r to update, respectively, for (19), the floating-point operations of U is (I 2 J + IJ 2 )R, and for (23), the floating-point operations of V is (IJ) 2 R. To sum up, the overall computational complexity of EWSP-NTF is O(IJKR + I 2 JR + IJ 2 R + I 2 J 2 R) at each iteration.
Second, the initialization of variable matrices is also a noticeable issue, which has a great impact on the performance of the method. In the experiments, We discuss two strategies: random initialization and VCA-FCLS initialization. The details of them can be seen in Section IV-A. Random initialization is used as a general strategy to initialize the abundance matrix E and the endmember matrix C in this article. That is, the values in the range [0, 1] are randomly selected to fill the entries of E and C, then the obtained results are taken as the initial matrices, respectively.
Third, stopping criterion is another issue that needs to be determined. We design two termination cases in our experiments: one is that the relative changes of abundance matrix and endmember matrix are less than the predefined error tolerance of 2 × 10 −4 for five times, and the other is that the maximum iteration number of 3000 is reached. When either of the above two cases is met, the optimization procedure stops.
Last but not the least, our proposed model EWSP-NTF is nonconvex, and meanwhile, it is hard to theoretically guarantee that the method is convergent with the adaptive weight factors. However, EWSP-NTF presents a stable convergence behavior, which is confirmed in subsequent experiments.
The unmixing performance use spectral angle distance (SAD), root mean square (RMSE), image reconstruction error (RE) as measured metrics. The SAD evaluates the dissimilarity between the real rth endmember signatureĉ r and its estimated signature c r . It is defined as The RMSE quantifies the difference between the real abundance mapÊ r of the rth endmember and its estimated map E r . It is defined as where N = I × J is the number of all pixels in HSI. Generally, there are no referenced endmembers and abundances in real data. The RE is used to measure the difference between the observed HSI Y and the reconstructed HSIŶ. It is defined as where K and N are the number of bands and the number of pixels in HSI, respectively. For all the above indices, the lower the value, the better the result.

A. Synthetic Data Experiments 1) Data Generation:
The synthetic data are generated by a similar way as in [39]. First, we select six pure spectral signatures from the United States Geological Survey digital spectral library as endmembers. Specifically, these are Carnallite, Ammoniojarosite, Almandine, Brucite, Axinite, and Chlonte, with 224 spectral bands ranging form 0.38 to 2.5 μm. The signatures of these selected endmembers are shown in Fig. 2.
Then, we generate the abundance as the following steps.
Step 1: Generate an image with size z 2 × z 2 pixels divided into z 2 blocks. Each block has z × z pixels.
Step 2: Full each pixel of a block with two kinds of randomly selected endmembers, and set their mixing ratios as θ and 1 − θ, respectively.
Step 4: Rescale the fractions of all endmembers in each pixel to satisfy the ANC and ASC constraints.
With endmembers and abundance maps constructed above, we can obtain a clean synthetic HSI by LMM. In order to further evaluate the robustness to additive noise, zero-mean white Gaussian noise is added to the clean HSI. Signal-to-noise ratio (SNR), the measure of noise level, is defined as follows: where y and e represent the clean data and noise data of a pixel, respectively. E[·] denotes the expectation operator. A series of generated noisy images with SNR will be used in later experiments.
2) Parameter Setting: In the proposed algorithm, there are four parameters δ, μ, λ, and ρ to be determined. In order to reflect the effect of each parameter on the model, we analyze the change of the parameter in a variation range by fixing other parameters for the clean synthetic data (where z = 8 and θ = 0.8).
First, we set δ ∈ {0, 0.2, 0.4, 0.6, 0.8, 1, 3, 5, 8, 10}, which gives the impact of ASC on the unmixing performance. Fig. 3(a)  shows the index results of SAD and RMSE when δ changes. From the figure, it can be seen that δ has a stable result between 0.6 and 1. We choose the optimal value δ = 0.8 in the following experiments.
Then, we discuss the balance parameter μ within the range of {20, 30, 40, 50, 60, 80, 100, 150}, which controls the convergence rate of our proposed unmixing algorithm. Fig. 3(b) shows the index results of SAD and RMSE when μ changes. According to the figure, both SAD and RMSE obtain the best results when μ = 50, thus this value is set in the following experiments.
Finally, as λ and ρ are the mutually influential parameters in our proposed model, which represent the lowrank and sparse constraints of the abundance map on the unmixing performance, respectively, we use a grid method with λ ∈ {0.001, 0.01, 0.1, 1, 50, 100, 250, 500} and  ρ ∈ {0.001, 0.01, 0.1, 1, 50, 100, 250, 500} to jointly search for their optimal values. Fig. 4(a) and (b) displays the changes of SAD and RMSE values, respectively. In these two figures, the change of ρ is more gentle than that of λ, and both of them have slight change between 0.001 and 1. For the following experiments, we set the coordinates of the concave point λ = 1 and ρ = 50 illustrated in the figures as their parameter values, which contribute to the best unmixing result.
3) Scene Variation: According to the synthesis method, there are three factors that can be debugged to generate the experimental simulated data, which are noise level SNR, image size z, and endmembers' mixing ratio θ. In addition, the number of endmembers is usually unknown in the actual scene. Thus, we discuss these four variables within the appropriate scopes in the following, and test the performance of the seven unmixing methods. The parameters of all the used methods are adjusted to their optimal values for experiments.  Fig. 5(b) shows the index results of all methods when z changes. This figure reflects the phenomenon that more pixels contain more abundant information, which is more conducive to unmixing. It is worth noting that our proposed method gets the best results in both SAD and RMSE of all the set image sizes.
Then, we evaluate the performance of different methods under different mixing levels. The mixing level is represented by ratio θ, where a smaller θ implies a higher mixing level. We construct the synthetic data by fixing SNR = 25 dB, z = 8, and θ is chosen from {0.5, 0.6, 0.7, 0.8, 0.9}. Fig. 5(c) presents the experimental results of all methods varying with the mixing level. As shown in this figure, all methods perform better when θ increases. Our proposed method has better results at higher θ values compared with other methods, which indicates that its advantages are enhanced in highly correlated data.
Finally, we analyze the impact of different methods under different estimated number of endmembers. The data of z = 8, θ = 0.8, SNR = 25 dB, and the true number of endmembers R = 6 are tested. We select the estimated number of endmembers in range of {4, 5, 6, 7, 8} for experiments. Fig. 5(d) shows the SAD and RMSE results of all methods adopting different estimates. It can be seen that the number of endmembers gets the best result at the exact value. The misestimation of the number of endmembers affects the unmixing performance of all methods.
However, our method can still keep the advantage in a small range of the misestimation.

4) Unmixing Performance:
We select a set of unmixing results for the synthetic data with the size of 64 × 64 × 224, θ = 0.8, and SNR = 25 dB for demonstration. Table I gives the numerical index results, in which the bold value corresponds to the best result of each endmember. Compared with these methods, our method performs better for endmembers #2, #3, #5, and #6 in terms of SAD. Besides, the endmembers except endmembers #2 and #4 also get the best RMSE. The mean values of SAD and RMSE indicate that our method is superior to other methods on the overall performance. Figs. 6 and 7 visually show the advantages of the endmember curves and abundance maps restored by our method.

5) Initialization Option:
We test the unmixing methods under two different initialization strategies: random initialization and VCA-FCLS initialization for the synthetic data with z = 8, θ = 0.8, and SNR = 25 dB. All the methods have 20 trials at the same initial value. The means and standard deviations of SAD and RMSE are shown in Fig. 8. This figure shows that all the methods have the better numerical results by VCA-FCLS strategy. And our method achieves the best performance on both initialization strategies. The advantage of our method in adopting random initialization is more obvious, which indicates that it is robust in initial selection. 6) Ablation Analysis: Here, we design two ablation experiments to discuss the role of regular terms and adaptive parameters in our model. First, the ablation experiment of regularizers is performed by ignoring the sum-to-one constraint, the low-rank constraint, and the sparse constraint in our proposed model, respectively. Corresponding to the following three loss functions:    We test them and the original loss function (6) for the synthetic data with z = 8, θ = 0.8, and SNR = 25 dB. Table II gives the SAD and RMSE values of 20 trials. It can be seen that the proposed method with all constraints has the best results. Figs. 9 and 10 show the endmember curves and abundance maps for endmember #3 obtained by the four models, respectively.  Comparing them with the references, we can draw the conclusion that each regularizer in our model has a positive role on the unmixing performance.
Then, we discuss the effect of the low-rank and sparse adaptive parameters α r and β r in our model (6). Fig. 11 shows the changes of adaptive parameters during the iteration for the synthetic data. It can be seen that both the two adaptive parameters of each abundance map tend to different stable values as iteration increases. This behavior confirms that each abundance map has different properties and it is necessary to treat it differently. Four cases are considered in the ablation experiment of adaptive parameters: with no adaptive parameters (when α r = 1, β r = 1), with sparse adaptive parameter (when α r = 1), with low-rank adaptive parameter (when β r = 1), and with both low-rank and sparse adaptive parameters (the original model). We test these cases for the 20 groups synthetic data with z = 8, θ = 0.8, and SNR = 25 dB. Table III gives the mean and standard deviation of SAD and RMSE. The results show that it is advantageous to have adaptive parameters.

B. Real-World Data Experiments
Four real-world datasets are used to test the performance of the proposed method: Urban, Jasper Ridge, Samson, and For visual comparison with other methods, we follow the same approach in [39] to get the referenced endmembers and abundances: the strategy of initialization and parameter debugging is as well as that adopted in the synthetic data experiments. 1) Urban Data: Urban dataset is one of the most widely used data in HU study, which is acquired by the hyperspectral digital imagery collection experiment (HYDICE) sensor over an urban area. This image has 307 × 307 pixels with 210 spectral bands, whose wavelength ranges from 0.4 to 2.5 μm. Affected by the imaging environment, we remove some low-SNR bands and water-vapor bands (1-4, 76, 87, 101-111, 136-153, and 198-210) from the original data, leaving 162 bands for experiments. Besides, we assume that Urban data are consisted of four endmembers: asphalt, grass, tree, and roof. Table IV presents the SAD and RMSE results of the seven unmixing methods. From Table IV, we see that MRS-NMF and TV-RSNMF obtain the best SAD for the endmember asphalt and the endmember grass, respectively. MV-NTF have the best RMSE for the endmember grass. And our proposed method achieves the best results in other endmembers, and the mean value of SAD and RMSE is far superior to other methods. Figs. 12 and 13 show the results of endmembers and abundances obtained by unmixing the Urban data. Fig. 12 displays that TV-RSNMF is closest to the reference curve for the endmember grass, and our method performs well in the remaining endmembers. According to Fig. 13, all methods can restore the general information compared with the referenced abundance maps. However, our method is superior in details since more properties of abundances are considered.
2) Jasper Ridge Data: Jasper Ridge dataset is another commonly used data for HSI processing, which is collected by the airborne visible/infrared imaging spectrometer sensor over the Jasper Ridge Natural Reserve in California, USA. The original data contain 512 × 614 pixels with 224 bands ranging from 0.38 to 2.5 μm. In order to ensure the efficiency of the experiments, we remove the useless bands (1-3, 108-112, 154-166, and 220-224) and intercept the data of 100 × 100 pixels with 198 bands for use. Four endmembers are assumed to exist in this subimage: tree, water, soil, and road.   see that MRS-NMF gets the best SAD for the endmember tree, and TV-MV-NTF achieves the best RMSE for the endmember road. Apart from these, the best results of other endmembers are obtained by our proposed method. The mean value of SAD and RMSE indicates that our method is better than other methods on overall unmixing performance. Figs. 14 and 15 visually show the difference between estimated and referenced endmembers and abundance maps for Jasper Ridge data, respectively. Fig. 14 displays that the endmember curves of MRS-NMF recovered well for endmembers tree, water, and soil, but the curve for the endmember road is quite different from the referenced curve. The recovery effect of our method for four endmembers is relatively stable. According to Fig. 15, the abundance maps restored by our proposed method are more consistent with the references, especially for endmembers tree, water, and soil.
3) Samson Data: Samson hyperspectral dataset is collected by Florida Environment Research Institute using Samson sensor. The full image contains 952 × 952 pixels with 156 bands, and its wavelength ranges from 0.401 to 0.889 μm. In our experiments, there is a part of Samson data with the size of 95 × 95 pixels for unmixing. This subimage is set to include three endmembers: soil, tree, and water. Table VI gives the SAD and RMSE values of the seven unmixing methods for Samson data. It can be seen that our    Fig. 16 gives the spectral curves that are referenced and restored by unmixing methods. It can be found that the curves obtained by our method have the smallest difference with the references. From the restored and corresponding referenced abundance maps shown in Fig. 17, our method gets the closest result to the references. It illustrates that the low-rank and sparse constraints added by our method to the abundance maps are effective.
Moreover, the RE is considered to compare the difference between the real ground data Y andŶ = R r=1Ê r •ĉ r reconstructed by unmixing components here. The index does not require the referenced endmembers and abundance maps in real data.   endmembers considered in this scene are grass, tree, road, roof, water, and trail.
Since the deep learning-based unmixing methods widely use this data for test, we also compare the performance of our method and three deep learning-based methods UnDIP [50], CyCUNet [49], MiSiCNet [51] here. Figs. 18 and 19 show the spectral curves and abundance maps obtained by different unmixing methods for endmembers tree and roof, respectively. It can be seen that UnDIP obtains the most similar unmixing result to the reference. Although our method lacks the specific training of data, it is still competitive compared with these deep learning-based methods and is superior to the same type of traditional unmixing methods for this data.

C. Discussion
Not only does EWSP-NTF improve the unmixing performance, but its novelty lying in the combination of the equivalence with NMF also opens up a new algorithm for the MV-NTF-based methods. Therefore, in this section, we record the running time of all methods to further compare their computational efficiency. And the stability of the proposed method under different experimental data is also demonstrated.
1) Computing Time: In order to tolerate the differences among methods, the stopping conditions of all methods keep their original optimal strategies in our experiments. Thus, for the fairness of time comparison, we record the per iteration time of methods under the same computer equipment. Specifically, Table VIII displays the time cost per iteration of different methods for the experimental datasets. As seen in Table VIII, NMF takes the shortest time because it has no constraints. Although our proposed method adds the abundance constraints to the basic MV-NTF framework, even the low-rank constraint needs the SVD, it still has a great speed-up among these MV-NTF-based methods. This time advantage can be explained by comparing their algorithm complexity. The simplification of our method is to avoid the low-rank matrix decomposition, which can be done without the rank of the abundance map.
2) Convergence Analysis: For further convincing the stability of our proposed method, we numerically analyze its convergence by recording the objective function values during the iteration both for the synthetic data and real-world data. Fig. 20(a) depicts the convergence curve for the synthetic data with z = 8, θ = 0.8, and SNR = 25 dB. And Fig. 20(b) shows the convergence curve for Samson data. It can be observed from the figures that their objective function values are monotonously decreasing and gradually becoming stable as the iteration number increases. Besides, other data have the similar behavior in the experiment, so they are omitted here.

V. CONCLUSION
In this article, we propose a new unmixing method named EWSP-NTF for the spectral unmixing problem. This method not only retains the lossless third-order tensor representation of HSI but also considers the characteristics of each abundance map, which are constrained by adaptive low-rank and sparse regular terms simultaneously. Compared with the existing methods, EWSP-NTF distinguishes to treat each abundance map, which is more in line with the actual scene. Furthermore, abandoning the low-rank matrix decomposition of MV-NTF but considering the abundance map as a whole also allows it to simplify the tensor operation by linking with NMF. Experiments on synthetic data and real-world data confirm the effectiveness of our proposed method. EWSP-NTF achieves the better unmixing performance and has advantages in time consumption against the same type of tensor-based unmixing methods. Besides, our method has a competitive performance with the deep learning-based unmixing methods. In the future, we will explore the potential properties of endmembers and abundances to further improve the unmixing performance. Moreover, we will study the unmixing application of NTF in nonlinear mixing mode.