Low-Rank Tensor Modeling for Hyperspectral Unmixing Accounting for Spectral Variability
Tales Imbiriba, Ricardo Augusto Borsoi, José Carlos Moreira Bermudez
AACCEPTED FOR PUBLICATION IN THE IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING – OCTOBER 25, 2019. 1
Low-Rank Tensor Modeling for HyperspectralUnmixing Accounting for Spectral Variability
Tales Imbiriba,
Member, IEEE , Ricardo Augusto Borsoi,
Student Member, IEEE ,Jos´e Carlos Moreira Bermudez,
Senior Member, IEEE
Abstract —Traditional hyperspectral unmixing methods neglectthe underlying variability of spectral signatures often observed intypical hyperspectral images (HI), propagating these missmodel-ing errors throughout the whole unmixing process. Attempts tomodel material spectra as members of sets or as random variablestend to lead to severely ill-posed unmixing problems. Althoughparametric models have been proposed to overcome this draw-back by handling endmember variability through generalizationsof the mixing model, the success of these techniques dependon employing appropriate regularization strategies. Moreover,the existing approaches fail to adequately explore the naturalmultidimensinal representation of HIs. Recently, tensor-basedstrategies considered low-rank decompositions of hyperspectralimages as an alternative to impose low-dimensional structures onthe solutions of standard and multitemporal unmixing problems .These strategies, however, present two main drawbacks: 1) theyconfine the solutions to low-rank tensors, which often cannotrepresent the complexity of real-world scenarios; and 2) they lackguarantees that endmembers and abundances will be correctlyfactorized in their respective tensors. In this work, we propose amore flexible approach, called ULTRA-V, that imposes low-rankstructures through regularizations whose strictness is controlledby scalar parameters. Simulations attest the superior accuracyof the method when compared with state-of-the-art unmixingalgorithms that account for spectral variability.
Index Terms —Hyperspectral data, endmember variability, ten-sor decomposition, low-rank, ULTRA, ULTRA-V.
I. I
NTRODUCTION
Hyperspectral imaging has attracted formidable interestfrom the scientific community in the past two decades, andhyperspectral images (HI) have been explored in a vast, andincreasing, number of applications in different fields [1]. Thelimited spatial resolution of hyperspectral devices often mixesthe spectral contributions of different pure materials, termed endmembers (EM), in the scene [2]. This phenomenon ismore explicit in remote sense applications, due to the distancebetween airborne and spaceborne sensors and the target scene.
This work has been supported by the National Council for Scien-tific and Technological Development (CNPq) under grants 304250/2017-1,409044/2018-0, 141271/2017-5 and 204991/2018-8, and by the BrazilianEducation Ministry (CAPES) under grant PNPD/1811213.T. Imbiriba was with the Department of Electrical Engineering, FederalUniversity of Santa Catarina (DEE–UFSC), Florian´opolis, SC, Brazil, andis with the ECE department of the Northeastern University, Boston, MA,USA. e-mail: [email protected]. R.A. Borsoi is with the DEE–UFSC,Florian´opolis, SC, Brazil, and with the Lagrange Laboratory, Universit´e Cˆoted’Azur, Nice, France. e-mail: [email protected]. J.C.M. Bermudez is withthe DEE–UFSC, Florian´opolis, SC, Brazil, and with the Graduate Programon Electronic Engineering and Computing, Catholic University of Pelotas(UCPel) Pelotas, Brazil. e-mail: [email protected] received Month day, year; revised Month day, year.
The mixing process must be well understood for the vital infor-mation relating the pure materials and their distribution in thescene to be accurately unveiled. Hyperspectral unmixing (HU)aims at solving this problem by decomposing the hyperspectralimage into a collection of endmembers and their fractional abundances [3].Different mixing models have been employed to explainthe interaction between light and the endmembers [1]-[4].The simplest and most widely used model is the LinearMixing Model (LMM) [2], which assumes that the observedreflectance vector ( i.e. a pixel) can be modeled as a convexcombination of the spectral signatures of each endmemberpresent in the scene. Convexity imposes positivity and sum-to-one constraints on the linear combination coefficients.Hence, they represent the fractional abundances with whichthe endmembers contribute to the scene. Though the simplicityof the LMM leads to fast and reliable unmixing strategiesin some situations, it turns out to be simplistic to explainthe mixing process in many practical applications. Hence,several approaches have been proposed in the literature toaccount for nonlinear mixing effects [4]–[6] and endmembervariability [7]–[9] often present in practical scenes.A myriad of factors can induce endmember variability,including environmental, illumination, atmospheric and tem-poral changes [7]. If not properly considered, such variabilitycan result in significant estimation errors being propagatedthroughout the unmixing process [10]. Most of the methodsproposed so far to deal with spectral variability can be classi-fied in three major groups: endmembers as sets, endmembersas statistical distributions and, more recently, methods thatincorporate the variability in the mixing model, often usingphysically motivated concepts [8]. The method proposed inthis work belongs to the third group. Parametric models havereceived considerable interest since, unlike the other twoapproaches, they require neither spectral libraries to be know apriori nor strong hypothesis about the endmembers statisticaldistribution.Recently, [10], [11] and [12] introduced variations of theLMM to cope with the spectral variability. Unmixing usingthese models lead to ill-posed problems that were solved byusing a combination of different regularizations terms andvariable splitting optimization strategies. The model PerturbedLMM (PLMM) in [10] augmented the endmember matrix withan additive perturbation matrix that needs to be estimatedjointly with the abundances. Although the additive perturbationcan model arbitrary endmember variations, it is not physicallymotivated, and the excessive amount of degrees of freedom a r X i v : . [ c s . C V ] O c t ACCEPTED FOR PUBLICATION IN THE IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING – OCTOBER 25, 2019. makes the problem even harder to solve. The Extended LMM(ELMM) proposed in [11] introduces one new multiplicativeterm for each endmember, and can efficiently model changesin the observed reflectance due to illumination effects [11].This model has a clear physical motivation, but its modelingcapability is limited. The Generalized Linear Mixing Model(GLMM) proposed in [12] generalizes the ELMM to accountfor variability in all regions of the measured spectrum. TheGLMM is physically motivated and capable of modeling arbi-trary variability, resulting in improved accuracy at the expenseof a small increase in the computational complexity, whencompared to the ELMM. Other works attempt to capture thesecomplex spectral variations indirectly by means of additiveresidual terms [13], [14]. Although avoiding the interactionsbetween the abundance fractions and the endmember signa-tures, these strategies usually do not estimate the EM spectrafor each image pixel.The above mentioned methods resort to different strategiesto regularize the ill-posed optimization problem leading to theestimation of abundances and endmembers. The regularizationis achieved by introducing into the unmixing problem addi-tional information based on common knowledge about thelow-dimensionality of structures embedded in hyperspectralimages.Possible ways to recover lower-dimensional structures fromnoisy and corrupted data include the imposition of low-rankmatrix constraints on the estimation process [15], or the low-rank decomposition of the observed data [16], [17]. The factsthat HIs are naturally represented and treated as tensors, andthat low-rank decompositions of higher-order ( ą
2) tensorstend to capture homogeneities within the tensor structuremake such strategies even more attractive for HU. Low-ranktensor models have been successfully employed in varioustasks involving HIs, such as recovery of missing pixels [18],anomaly detection [19], classification [20], compression [21],dimensionality reduction [22] and analysis of multi-angleimages [23]. More recently, [24] and [23] considered low-ranktensor decompositions applied to standard and multitemporalHU, respectively.In [24] the HI is treated as three-dimensional tensor, andspatial regularity is enforced through a nonnegative tensorfactorization (NTF) strategy that imposes a low-rank tensorstructure. In [23], nonnegative canonical polyadic decompo-sition were used to unmix multitemporal HIs represented asthree-dimensional tensors built by stacking multiple temporalmatricized HIs. Though a low-rank tensor representation maynaturally describe the regularity of HIs and abundance maps,the forceful introduction of stringent rank constraints mayprevent an adequate representation of fast varying structuresthat are important for accurate unmixing. Another limitationof the approaches proposed in [24] and [23] is the lack ofguarantee that endmembers and abundances will be correctlyfactorized into their respective tensors. In [25], we proposeda new low-rank HU method called
Unmixing with Low-rankTensor Regularization Algorithm (ULTRA), which accountsfor highly correlated endmembers. The HU problem wasformulated using tensors and a low-rank abundance tensorregularization term was introduced. Differently, from the strict tensor decomposition considered in [23], [24], ULTRAallowed important flexibility to the rank of the estimatedabundance tensor to adequately represent fine scale structureand details that lie beyond a low-rank structure, but withoutcompromising the regularity of the solution.In this work we extend the strategy proposed in [25] toaccount for the important effect of endmember variability aswell as a novel method to estimate the sufficient rank of atensor for accurately solving the HI unmixing problem. Themain novel contributions of this paper are:a) We extend the strategy proposed in [25] by imposing anew low-rank regularization on the four-dimensional end-member tensor, which contains one endmember matrixfor each pixel, to account for endmember variability. Thenew cost function results in an iterative algorithm, named
Unmixing with Low-rank Tensor Regularization Algo-rithm accounting for endmember Variability (ULTRA-V).At each iteration, ULTRA-V updates the estimations ofthe abundance and endmember tensors as well as theirlow-rank approximations.b) We propose a novel non-trivial strategy to determine thesmallest rank representation that contains most of thevariation of multilinear singular values [26].Simulation results using synthetic and real data illustrate theperformance improvement obtained using ULTRA-V whencompared to competing methods, as well as its competitivecomputational complexity for relatively small images.The paper is organized as follows. Section II briefly reviewsimportant background on linear mixing models and definitionsand notation used for tensors. Section III presents the proposedsolution and the strategy to estimate tensor ranks. Section IVpresents the simulation results and comparisons. Finally, Sec-tion V presents the conclusions.II. B
ACKGROUND AND NOTATION
A. Extended Linear Mixing Models
The Linear Mixing Model (LMM) [2] assumes that a givenpixel r n “ r r n, , . . . , r n,L s J , with L bands, is represented as r n “ M α n ` e n , subject to J α n “ and α n ľ (1)where M “ r m , . . . , m R s is an L ˆ R matrix whosecolumns are the R endmembers m i “ r m i, , . . . , m i,L s J , α n “ r α n, , . . . , α n,R s J is the abundance vector, e „ N p , σ n I L q is an additive white Gaussian noise (WGN), I L isthe L ˆ L identity matrix, and ľ is the entrywise ě operator.The LMM assumes that the pure material endmembers arefixed for all pixels r n , n “ , . . . , N , in the HI. This assump-tion can jeopardize the accuracy of estimated abundances inmany circumstances due to the spectral variability existing ina typical scene.Different extensions of the LMM have been recently pro-posed to mitigate this limitation. These models employ adifferent endmember matrix for each pixel, and are particularcases of the model r n “ M n α n ` e n (2) MBIRIBA, BORSOIM, BERMUDEZ 3 where M n is the endmember matrix for the n -th pixel.Different parametric models propose different forms for M n to account for spectral variability. These include additiveperturbations over a mean matrix in the PLMM [10], mul-tiplicative factors applied individually to each endmember inthe ELMM [11] or to each band in the GLMM [12]. Moreover,spatial regularization of the multiplicative scaling factors in theELMM and GLMM help to further mitigate the ill-posednessof the problem. Note that, although some works proposedto handle complex spectral variations indirectly by means ofadditive residual terms [13], [14], these strategies usually donot estimate the EM spectra for each image pixel and, like theother models, also require carefully designed regularizationstrategies.Other approaches pursue different ways to improve theconditioning of the inverse problem, employing for instancemultiscale regularization on the abundance maps [27], or usingadditional information in the form of spectral libraries knowna priori [28], [29] or extracted from the observed HI [30].All these models, however, fail to exploit the high dimen-sional structure of the problem, which naturally suggests therepresentation of the HI, abundance maps, and endmembermatrices for all pixels as higher-order tensors. In this work,instead of introducing a rigid parametric model for the end-members, we employ a more general tensor model, using awell-devised low-rank constraint to introduce regularity to theestimated endmember tensor. B. Notation
An order- P tensor T P R N ˆ¨¨¨ˆ N P ( P ą ) is an N ˆ ¨ ¨ ¨ ˆ N P array with elements indexed by T n ,n ,...,n P .The P dimensions of a tensor are called modes . A mode- (cid:96) fiber of tensor T is the one-dimensional subset of T obtained by fixing all but the (cid:96) -th dimension, and is indexedby T n ,...,n (cid:96) ´ , : ,n (cid:96) ` ,...,n P . A slab or slice of tensor T isa two-dimensional subset of T obtained by fixing all buttwo of its modes. An HI is often conceived as a threedimensional data cube, and can be naturally represented byan order-3 tensor R P R N ˆ N ˆ L , containing N ˆ N pixelsrepresented by the tensor fibers R n ,n , : P R L . Analogously,the abundances can also be collected in an order-3 tensor A P R N ˆ N ˆ R . Thus, given a pixel R n ,n , : , the respectiveabundance vector α n ,n is represented by the mode-3 fiber A n ,n , : . Similarly, the endmember matrices for each pixelcan be represented as an order-4 tensor M P R N ˆ N ˆ L ˆ R ,where M n ,n , : , : “ M n ,n . We now review some operationsof multilinear algebra (the algebra of tensors) that will be usedin the following sections (more details can be found in [31]). C. Tensor product definitions
Definition 1.
Outer product:
The outer product between vec-tors b p q P R N , b p q P R N , . . . , b p P q P R N P is defined as theorder- P tensor T “ b p q ˝ b p q ˝ ¨ ¨ ¨ ˝ b p P q P R N ˆ N ˆ¨¨¨ˆ N P ,where T n ,n ,...,n P “ b p q n b p q n ¨ ¨ ¨ b p P q n P and b p i q n i is the n i -thposition of b p i q . It generalizes the outer product between twovectors. Fig. 1: Polyadic decomposition of a three-dimensional tensor,written as both outer products and mode ´ n products. Definition 2.
Mode- k product: The mode- k product, denoted U “ T ˆ k B , of a tensor T P R N ˆ¨¨¨ˆ N k ˆ¨¨¨ˆ N P and amatrix B P R M k ˆ N k is evaluated such that each mode- k fiberof T is multiplied by matrix B , yielding U n ,n ,...,m k ,...,n P “ ř N k i “ T ...,n n ´ ,i,n n ` ,... B m k ,i . Definition 3.
Multilinear product:
The full multilinear prod-uct, denoted by T ; B p q , B p q , . . . , B p P q , consists of thesuccessive application of mode- k products between T and ma-trices B p i q , represented as T ˆ B p q ˆ B p q ˆ . . . ˆ P B p P q . Definition 4.
Mode-(M,1) contracted product:
The con-tracted mode-M product, denoted by U “ T ˆ M b , is aproduct between a tensor T and a vector b in mode-M,where the resulting singleton dimension is removed, given by U ...,n n ´ ,n n ` ,... “ ř N n i “ T ...,n n ´ ,i,n n ` ,... b i .D. The Canonical Polyadic Decomposition An order- P rank-1 tensor is obtained as the outer productof P vectors. The rank of an order- P tensor T is defined asthe minimum number of order- P rank-1 tensors that must beadded to obtain T [17]. Thus, any tensor T P R N ˆ N ˆ¨¨¨ˆ N P with rank p T q “ K can be decomposed as a linear combina-tion of at least K outer products of P rank-1 tensors. This so-called polyadic decomposition is illustrated in Figure 1. Whenthis decomposition involves exactly K terms, it is called thecanonical polyadic decomposition (CPD) [31] of a rank- K tensor T , and is given by T “ K ÿ i “ λ i b p q i ˝ b p q i ˝ ¨ ¨ ¨ ˝ b p P q i . (3)It has been shown that this decomposition is essentiallyunique under mild conditions [17]. The CPD can be writtenalternatively using mode- k products as T “ D λ ˆ B p q ˆ B p q ¨ ¨ ¨ ˆ P B p P q (4)or using the full multilinear product as T “ D λ ; B p q , B p q , . . . , B p P q (5)where D λ “ Diag P ` λ , . . . , λ K ˘ is the P -dimensional di-agonal tensor and B p p q “ r b p p q , . . . , b p p q K s , for p “ , . . . , P .Given a tensor T P R N ˆ N ˆ¨¨¨ˆ N P , the CPD can be obtainedas the solution to the following optimization problem [17] min D λ , B p q ,..., B p K q ››› T ´ K ÿ i “ λ i b p q i ˝ ¨ ¨ ¨ ˝ b p P q i ››› F . (6)A widely used strategy to compute an approximate solutionto (6) is to use an alternating least-squares technique [17], ACCEPTED FOR PUBLICATION IN THE IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING – OCTOBER 25, 2019. which optimizes the cost function with respect to one termat a time, while keeping the others fixed, until convergence.Although optimization problem (6) is generally non-convex,its solution is unique under relatively mild conditions, whichis an important advantage of tensor-based methods [17].
E. Tensor rank bounds
Finding the rank of an arbitrary tensor T is NP-hard [32].In [17], upper and lower bounds on tensor ranks are presentedfor arbitrary tensors. Let T be an order-3 tensor and R ” dim span t T : ,j,k u @ j,k R ” dim span t T i, : ,k u @ i,k (7) R ” dim span t T i,j, : u @ i,j be the mode-1 (column), mode-2 (row) and mode-3 (fiber)ranks, respectively, of T . Thus, the K ” rank p T q which isable to represent an arbitrary tensor is limited in the interval max p R , R , R q ď K ď min p R R , R R , R R q . (8)The reader can note that the bounds presented above oftenlead to very large tensor ranks. In many practical applications,however, the “useful signal” rank is often much less than theactual tensor rank [17]. Hence, when low-rank decompositionsare employed to extract low-dimensional structures from thesignal, the ranks that lead to meaningful results are usuallymuch smaller than max p R , R , R q . In Section III-F, wepropose a strategy to estimate the rank of tensor CPDs basedon the variation of the multilinear singular values of T .III. L OW - RANK U NMIXING P ROBLEM
An effective strategy to capture the low-dimensional struc-tures of HIs for solving the HU problem is to impose alow-rank structure to the abundance tensor [24]. The samestrategy can also be applied to the endmember tensor if oneconsiders the endmember variabilities to be small or highlycorrelated in low-dimensional structures within the HI. Thelow-rank property of HI tensors has been an important toolin the design of hyperspectral image completion [33] andrestoration algorithms [34], consisting in one of the main low-dimensional structures that are currently being considered inhyperspectral imaging applications. Thus, assuming that A hasa low-rank K A , and that M has a low-rank K M the globalcost functional for the unmixing problem can be written as J p A , M q “ N ÿ n “ N ÿ n “ } R n ,n , : ´ M n ,n , : , : A n ,n , : } F s. t. rank p M q “ K M , M ľ (9)rank p A q “ K A , A ľ , A ˆ R “ N ˆ N . Defining the HU problem as in (9) with fixed data independentranks K M and K A limits its flexibility to adequately repre-sent the desired abundance maps and endmember variability.Though fixing low ranks for A and M tends to capture themost significant part of the tensors energy [35], one may incurin a loss of fine and small scale details that may be relevant forspecific data. On the other hand, using large values for K A and K M makes the solution sensitive to noise, undermining thepurpose of regularization. Thus, an important issue is how toeffectively impose the low-rank constraint to achieve regularityin the solution without undermining its flexibility to adequatelymodel small variations and details.We propose to modify (9) by introducing new regularizationterms, controlled by two low-rank tensors Q P R N ˆ N ˆ R and P P R N ˆ N ˆ L ˆ R , to impose non-strict constraints on K A and K M . Doing that, tensors Q and P work as a priori information, and the strictness of the low-rank constraint iscontrolled by two additional parameters λ A , λ M P R ` . Theproposed cost function is given by J p A , M , P , Q q “ N ÿ n “ N ÿ n “ } R n ,n , : ´ M n ,n , : , : A n ,n , : } F ` λ M } M ´ P } F ` λ A } A ´ Q } F s. t. M ľ , A ľ , A ˆ R “ N ˆ N (10)with rank p P q “ K M and rank p Q q “ K A . The optimizationproblem becomes p p A , x M , p P , p Q q “ arg min A , M , P , Q J p A , M , P , Q q . (11)To solve (11), we propose to find a local stationary pointby minimizing (10) iteratively with respect to each variable.The resulting algorithm is termed the Unmixing with Low-rank Tensor Regularization Algorithm accounting for spectralVariability (ULTRA-V), and is presented in Algorithm 1. Theintermediate steps are detailed in the following.
Algorithm 1:
Global algorithm for solving (10)
Input : R , λ M , λ A , A p q , and M p q . Output: p A and x M . K Q = estimateTensorRank( A p q ); K P = estimateTensorRank( M p q ); Set i “ ; while stopping criterion is not satisfied do i “ i ` ; P p i q “ arg min P J p A p i ´ q , M p i ´ q , P q ; Q p i q “ arg min P J p A p i ´ q , M p i ´ q , Q q ; M p i q “ arg min M J p A p i ´ q , M , P p i q , Q p i q q ; A p i q “ arg min A J p A , M p i q , P p i q , Q p i q q ; end return p A “ A p i q , x M “ M p i q ; A. Solving with respect to A To solve problem (11) with respect to A we use only theterms in (10) that depend on A , leading to the cost function J p A q “ N ÿ n “ N ÿ n “ } R n ,n , : ´ M n ,n , : , : A n ,n , : } F ` λ A } A ´ Q } F s. t. A ľ , A ˆ R “ N ˆ N (12) MBIRIBA, BORSOIM, BERMUDEZ 5 which results in a standard regularized fully constrained least-squares problem that can be solved efficiently.
B. Solving with respect to M Analogously to the previous section, to solve problem (11)with respect to M , we use only the terms in (10) that dependon M , leading to J p M q “ N ÿ n “ N ÿ n “ } R n ,n , : ´ M n ,n , : , : A n ,n , : } F ` λ M } M ´ P } F (13)s. t. M ľ . which results in a regularized nonnegative least-squares prob-lem. An approximate solution can be obtained ignoring thepositivity constraint over the endmember tensor and projectingthe least-squares result onto the positive orthant as [12] x M n ,n , : , : “ P ` ˆ´ R n ,n , : A J n ,n , : ` λ M P n ,n , : , : ¯´ A n ,n , : A J n ,n , : ` λ M I ¯ ´ ˙ (14)where P ` : R N ˆ N ˆ L Ñ R N ˆ N ˆ L ` is the projectionoperator that maps every negative element to zero. Althoughthis solution is approximate, it is significantly faster thandirectly solving (13) and the algorithm still demonstrated goodempirical convergence in our experiments. C. Solving with respect to P Rewriting the terms in (10) that depend on P leads to J p P q “ λ M } M ´ P } F . (15)Assuming that most of the energy of M lies in a low-rankstructure, we write the tensor P as a sum of a small number K P of rank-1 components, such that P “ K P ÿ i “ δ i x p q i ˝ x p q i ˝ x p q i ˝ x p q i . (16)This introduces a low-rank a priori condition on P whosestrictness can be controlled by the regularization constant λ P .Using (16) in (15) leads to the optimization problem ´ p ∆ , x X p q , x X p q , x X p q , x X p q ¯ “ (17) arg min ∆ , X p q , X p q , X p q , X p q λ M ››››› M ´ K P ÿ i “ δ i x p q i ˝ x p q i ˝ x p q i ˝ x p q i ››››› F where ∆ “ Diag p δ , . . . , δ K P q is a 4 dimensional diagonaltensor with ∆ i,i,i,i “ δ i . Problem (17) can be solved using analternating least-squares strategy [17].Finally, the solution p P is obtained from p ∆ , x X p q , x X p q , x X p q , and x X p q using the full multilinear product as p P “ ∆ ; x X p q , x X p q , x X p q , x X p q . (18) D. Solving with respect to Q Analogous to the previous section, the cost function to beoptimized for Q can be written as J p Q q “ λ A } A ´ Q } F . (19)Assuming that most of the energy of A lies in a low-rankstructure, we write tensor Q as a sum of a small number K Q of rank-1 components, such that Q “ K Q ÿ i “ ξ i z p q i ˝ z p q i ˝ z p q i . (20)This introduces a low-rank a priori condition on A , whichwill be more or less enforced depending on the regularizationconstant λ A . Using (20) in (19) leads to the optimizationproblem ´ p Ξ , p Z p q , p Z p q , p Z p q ¯ “ (21) arg min Ξ , Z p q , Z p q , Z p q λ A ››› A ´ K Q ÿ i “ ξ i z p q i ˝ z p q i ˝ z p q i ››› F where Ξ “ Diag ` ξ , . . . , ξ K Q ˘ is an order-3 diagonal tensorwith Ξ i,i,i “ ξ i . Problem (21) can be solved using analternating least-squares strategy [17]. Finally, the solution p Q is obtained from p Ξ , p Z p q , p Z p q and p Z p q using the fullmultilinear product as p Q “ Ξ ; p Z p q ; p Z p q ; p Z p q . (22) E. Computational complexity of Algorithm 1
The computational complexity of each iteration of Algo-rithm 1 can be measured as follows. The optimizations w.r.t. A and M both consist of regularized constrained LS problemswith N N R and N N LR variables respectively. Thus, theseproblems can be solved with a complexity of O ` p N N R q ˘ and O ` p N N LR q ˘ , respectively. The optimizations w.r.t.variables P and Q consist of CP decompositions of thesetensors with ranks K P and K Q , respectively. Considering analternating least squares (ALS) approach for the CPD, theseoptimization problems will have computational complexitiesof O ` K iter K P N N LR ˘ and O ` K iter K Q N N R ˘ , respec-tively, where K iter is the number of ALS iterations [36]. Thus,the overall complexity of the algorithm scales linearly withthe number of ALS iterations and with the tensor ranks, andcubically in the problem dimensions. When processing largedatasets, the extra complexity could be partially mitigated byapplying image segmentation or band selection [6], [37], [38]strategies. This analysis is beyond the scope of the presentwork and will be addressed in the future. F. Estimating tensor ranks
In Section II-E we have recalled important results relatingbounds for order-3 tensor ranks to the span of the matricizedversions of tensors. We have also noted, from our own expe-rience, that those bounds tend to indicate tensor ranks that arelarger than the rank associated with the information relevant
ACCEPTED FOR PUBLICATION IN THE IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING – OCTOBER 25, 2019. for HU. Our interest in HU is to model low-dimensionalstructures of the HI using low-rank tensors. At the same time,this low-rank representation should be rich enough to includeall dimensions of the original HI tensor that contain relevantinformation. Therefore, although the literature presents manyrank estimation strategies (see [39] and references therein), inthis work we exploit the rank bounds discussed in Section II-Eto approximate the “useful rank” of a tensor by the number ofthe largest singular values of their matricized versions requiredto represent most of the tensor energy.Let T i “ mat i p T q P R N i ˆp N ...N i ´ N i ` ,...N P q be the ma-tricization of an arbitrary tensor T P R N ˆ N ... ˆ N P obtainedby stacking all tensor fibers along the i -th tensor dimension.Let s i “ SVD p T i q be the set of singular values of T i , sorteddescending in value. Also, let d i “ diff p s i q be the vectorof first order differences of the elements of s i , such that, d p i q j “ s p i q j ´ s p i q j ` . Then, we define the i -th candidate forrank of T as the smallest index j such that | d p i q j | sufficientlysmall, namely, ˆ R i “ min j, s.t. , | d p i q j | ă ε (23)where ε is a parameter limiting the singular value variation.In all experiments reported here we used ε “ . . We haveexperimentally verified that the resulting abundance MSE hasvery low sensitivity to the choice of ε . Finally, we approximatethe rank of tensor T as K “ max t ˆ R , . . . , ˆ R P u . (24)For the experiments reported in this paper, we have useddefinition (24) to estimate K P and K Q in (16) and (20)from the abundance and endmember tensors estimated usingsimple unmixing strategies such as the scaled constrained leastsquares (SCLS) [11].IV. S IMULATIONS
In this section, the performance of the proposed method-ology is illustrated through simulations with both syntheticand real data. We compare the proposed ULTRA-V methodwith the the fully constrained least squares (FCLS), the scaledconstrained least squares (SCLS) [11], the PLMM [10], theELMM [11], and the GLMM [12]. To highlight the differ-ences between ULTRA-V and ULTRA [25], we also considerULTRA for simulations with synthetic data.To measure the accuracy of the unmixing methods weconsider the Mean Squared Error (MSE)MSE X “ N X } vec p X q ´ vec p p X q} (25)where vec p¨q is the vectorization operator X Ñ x , R a ˆ b ˆ c ÞÑ R abc , N X “ abc , and the Spectral Angle Mapper for the HISAM R “ N N ÿ n “ arccos ˆ r J n p r n } r n }} p r n } ˙ (26)and for the endmembers tensorSAM M “ N N ÿ n “ R ÿ k “ arccos ˜ m J k,n x m k,n } m k,n }} x m k,n } ¸ . (27) TABLE I: Simulation results using synthetic data. Data Cube 0 – DC0MSE A MSE M SAM M MSE R SAM R TimeFCLS 1.81 - - 16.11 5.92 0.42SCLS 0.68 175.79 6.19 59.38 5.42 0.38PLMM 0.76 - - 10.93 3.94 4.44Data Cube 3 – DC2FCLS 1.90 - - 2.03 13.69 0.30SCLS 0.71
All the algorithms were implemented in Matlab on a desktopcomputer equipped with an Intel Core I7 processor with4.2Ghz and 16Gb of RAM. In all cases, we used endmembersextracted using the VCA [40] either to build the referenceendmember matrix or to initialize the different methods, withthe number of endmembers R assumed to be known a priori.The abundance maps were initialized using the maps estimatedby the SCLS. A. Synthetic data
For a comprehensive comparison among the different meth-ods we created three synthetic datasets, namely Data Cube0 (DC0), Data Cube 1 (DC1) and Data Cube 2 (DC2),with 50 ˆ
50 pixels (DC0 and DC2) and 70x70 pixels (DC1).DC0 and DC1 were built using three 224-band endmembersextracted from the USGS Spectral Library [41], while DC2was built using three 16-band minerals often found in bodiesof the Solar System. For the three datasets, spatially correlatedabundance maps were used, as depicted in the first columnof Fig. 2. For DC0, we adopted the variability model usedin [11] (a multiplicative factor acting on each endmember).For DC1, we used the variability model according to thePLMM [10]. For DC2, we used the Hapke model [42] devisedto realistically represent the spectral variability introduceddue to changes in the illumination conditions caused by thetopography of the scene [11]. White Gaussian noise was addedto all datasets to yield a 30dB SNR.To select the optimal parameters for each algorithm,we performed grid searches for each dataset. We used
MBIRIBA, BORSOIM, BERMUDEZ 7
Fig. 2: Abundance maps of (top–down) DC0, DC1, and DC2for all tested algorithms. Abundance values represented bycolors ranging from blue ( α k “ ) to red ( α k “ ).parameter search ranges based on the ranges tested anddiscussed by the authors in the original publications. Forthe PLMM we used γ “ , since the authors fixed thisparameter in all simulations, and searched for α and β in the ranges r ´ , ´ , . , . , . , . , , s and r ´ , ´ , ´ , ´ s , respectively. For both ELMMand GLMM, we selected the parameters among thefollowing values: λ S , λ M P r . , . , , , , s , λ A P r ´ , ´ , . , . , . , , , and λ ψ , λ Ψ P r ´ , ´ , ´ s , while for the proposedULTRA-V we selected the parameters in theintervals λ A P r . , . , . , , , s and λ M P r . , . , . , . , . , s . For the ULTRA wesearched λ A in the same interval used for the ULTRA-V.The results are shown in Table I, were the best and secondbest results for each metric are marked in bold red andbold blue, respectively. ULTRA-V clearly outperformed thecompeting algorithms for all datasets in terms of MSE A . Forthe other metrics, the best results depended on the datasets.In terms of MSE M and SAM M , ULTRA-V yielded thebest results for DC0 and DC1. Finally, ULTRA-V resultsfor MSE R and SAM R were the second best for DC0 andthe best for DC2. The execution times, shown in the right-most columns of Table I, show that ULTRA-V requiredthe smallest execution time among the more sophisticatedalgorithms (PLMM, ELMM and GLMM) for DC0 and DC2,and comparable execution time for DC1. As expected, theULTRA method provided results that were often better thanthe FCLS but significantly worse than those obtained frommethods accounting for EM variability. This happens because ULTRA imposes a low-rank structure over the abundances butkeeps the EMs fixed for all pixels, what greatly limits thealgorithm capacity to adapt to EM variations along the image. A -3 M S E
68 0.5 M M S E M A
50 0100 M S E M A
500 0
Fig. 3: Parameters sensitivity to changes around the optimalvalues. Left: DC0, Middle: DC1, and Right: DC2.
1) Parameters sensitivity:
While we have proposed a strat-egy to determine rank values for tensors P and Q , the parame-ters λ M and λ A need to be selected by the user. We now studythe sensitivity of the ULTRA-V performance to variations ofthe parameters within the parameter search intervals presentedin the previous section. Figure 3 shows the values of MSE A resulting from unmixing the data using each combination ofthe parameter values. The sensitivity clearly tends to increasewhen values less than 1 are used for both parameters. Ourpractical experience indicates that good MSE A results can beobtained using λ M in r , s , and large values about 100 for λ A . Moreover, some insensitivity is verified for small changesin λ A about large values. Thus, searching λ A in r . , s with values spaced by decades as done for the examples inthe previous section seems reasonable.
2) Discussion:
A close look at Fig. 2 reveals the abun-dance maps estimated using ULTRA-V look noisier thanthose obtained using ELMM and GLMM. This is becausethe proposed approach does not impose the local smoothnessimposed by total variation (TV), but emphasizes abundanceregularity by enforcing a low-rank property. We note, however,that the spacial smoothness imposed by TV is not necessarilymandatory for a good abundance estimation, as can be verifiedfrom the results in Table I.
B. Real data
For the simulations with real data, we considered threedatasets, consisting of the Houston, Samson and Jasper Ridgeimages. All datasets were captured by the AVIRIS, whichoriginally has 224 spectral bands. For all images, the waterabsorption bands were removed resulting in 188 bands for theHouston image, 156 bands for the Samson image and 198bands for the Jasper Ridge image. The Houston data set isknown to have four predominant endmembers [11], [43]. TheSamson and Jasper Ridge images are known to have threeand four endmembers, respectively [44]. For all images theendmembers were extracted using the VCA [40]. Fig. 4 showsthe reconstructed abundance maps for all images and for alltested methods. The quantitative results are shown in Table II.Note that since the ground truth (correct) abundance values arenot available for these images, only the reconstruction errorMSE R has been used as a sort of quality verification.The last column in Fig. 4 shows that the proposed ULTRA-V method provided an accurate abundance estimation, clearly ACCEPTED FOR PUBLICATION IN THE IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING – OCTOBER 25, 2019.
FCLS V e g e t a t i o n SCLSU PLMM ELMM GLMM ULTRA-V M e t . R oo f s C o n c r e t e A s ph a l t FCLS W a t e r SCLSU PLMM ELMM GLMM ULTRA-V T r ee S o il FCLS T r ee s SCLSU PLMM ELMM GLMM ULTRA-V S o il W a t e r R oa d Fig. 4: Abundance maps of the Houston (upper panel), Samson(middle painel), and Jasper Ridge (bottom painel) data sets forall tested algorithms. Abundance values represented by colorsranging from blue ( α k “ ) to red ( α k “ ).outperforming the competing algorithms . In fact, for theConcrete and Metallic Roofs endmembers, the ULTRA-Vabundance map presents stronger Concrete and Metallic Roofscomponents in the stadium stands and stadium towers, respec-tively, when compared with the other methods equipped fordealing with spectral variability. The performance improve-ment provided by ULTRA-V is clearer for the Samson andJasper Ridge images. For instance, there is significantly lessconfusion between the Water, Tree and Soil endmembers inthe ULTRA-V results for the Samson image when compared tothose of the PLMM, ELMM, and GLMM methods. Similarly,the ULTRA-V reconstructed abundance maps of the JasperRidge image show a much stronger Water component in theriver and less confusion between the Tree, Soil and Roadendmembers.The objective metrics presented in Table II indicate thatULTRA-V yields competitive reconstruction errors in termsof MSE. These results, however, should be interpreted withproper care, as the connection of reconstruction error andabundance estimation is not straightforward. The differences between the ELMM and ULTRA-V results are lesssignificant for the Houston image.
The execution times in Table II indicate that, as discussedin Section III-E, ULTRA-V did not scale well with thelarger image sizes and higher number of endmembers, whichdirectly impacted the CPD stage of ULTRA-V. Moreover,the more complex images resulted in higher rank estimatesusing the strategy discussed in Section III-F. This indicatesthat there is still room for improving the proposed methodby either providing a segmentation strategy or using fasterCPD methods. This, however, is an open problem that will beaddressed in future works.To assess the estimated endmember variabilities, we an-alyzed the results for the Samson data set. We consideredtwo approaches. The first approach consisted in averagingthe projection of the estimated endmembers on the threeeigenvectors associated to the three largest eigenvalues foreach endmember. The results are shown in Fig. 7. Theseplots illustrate the endmember variances for each pixel, withred implying a large variance and blue a small variance.The second approach consisted in directly comparing theendmembers estimated with ULTRA-V and VCA. The resultsare shown in Fig. 8. These figures illustrate the ability of theproposed method to characterize the spectral variability whileenforcing a spatial structure for the estimated endmembers.To illustrate the role of the low-rank tensors P and Q , wecompare them to the abundances A and endmembers M inFigures 5 and 6, for the Jasper Ridge dataset. Figure 5 showsthe estimated abundances A and their low-rank counterpart Q .One can verify that Q (bottom row) has a very coarse spatialdistribution when compared with A (top row). This showsthat imposing the low-rank structure through a regularizationconstraint gives the resulting abundances enough flexibilityto model fine-scale spatial details while maintaining most ofits spatial distribution. Figure 6 leads to similar conclusionsfor the endmembers. One can note that the low-rank tensor P has a coarser structure when compared with the estimatedendmembers M . This distinction can be seen very clearly forthe Water and Road endmembers. Q Trees A Soil Water Road
Fig. 5: Comparison of tensors A and Q after ULTRA-Vconvergence for the Jasper Ridge data set.V. C ONCLUSIONS
In this paper, we proposed a new low-rank regulariza-tion strategy for introducing low-dimensional spatial-spectralstructure into the abundance and endmember tensors forhyperspectral unmixing considering spectral variability. Theresulting iterative algorithm, called ULTRA-V, imposes low-rank structure by means of regularizations that force most of
MBIRIBA, BORSOIM, BERMUDEZ 9
Fig. 6: Comparison of tensors M and P after ULTRA-Vconvergence for the Jasper Ridge data set. Water Tree Soil
Fig. 7: Average of the ULTRA-V endmembers tensor projec-tion over the 3 principal components for the Samson data set.TABLE II: Real Data.
Algorithm Houston Data Samson Jasper RidgeMSE R Time MSE R Time MSE R TimeFCLS 0.2283 1.90 0.0177 1.38 0.3567 1.59SCLS 0.0037 2.04 0.0041 1.29 0.0271 1.79PLMM 0.0190 454.90 0.0034 105.36 0.0257 72.86ELMM the energy of the estimated abundances and endmembers tolay within a low-dimensional structure. The proposed approachdoes not confine the estimated abundances and endmembersto a strict low-rank structure, which would not adequately ac-count for the complexity experienced in real-world scenarios.The proposed methodology includes also a strategy to estimatethe rank of the regularization tensors P and Q , leaving onlytwo parameters to be adjusted within a relatively reducedsearch space. Simulation results using both synthetic and realdata showed that the ULTRA-V can outperform state-of-the-art unmixing algorithms accounting for spectral variability.R EFERENCES[1] J. M. Bioucas-Dias, A. Plaza, G. Camps-Valls, P. Scheunders,N. Nasrabadi, and J. Chanussot, “Hyperspectral remote sensing dataanalysis and future challenges,”
IEEE Geoscience and Remote SensingMagazine , vol. 1, no. 2, pp. 6–36, 2013.
Fig. 8: VCA result (black) and ULTRA-V (gray) endmembersfor each pixel of the Samson data set. [2] N. Keshava and J. F. Mustard, “Spectral unmixing,”
IEEE SignalProcessing Magazine , vol. 19, no. 1, pp. 44–57, 2002.[3] 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 ofSelected Topics in Applied Earth Observations and Remote Sensing ,vol. 5, no. 2, pp. 354– 379, 2012.[4] N. Dobigeon, J.-Y. Tourneret, C. Richard, J. C. M. Bermudez,S. McLaughlin, and A. O. Hero, “Nonlinear unmixing of hyperspectralimages: Models and algorithms,”
IEEE Signal Processing Magazine ,vol. 31, no. 1, pp. 82–94, Jan 2014.[5] T. Imbiriba, J. C. M. Bermudez, C. Richard, and J.-Y. Tourneret,“Nonparametric detection of nonlinearly mixed pixels and endmemberestimation in hyperspectral images,”
IEEE Transactions on Image Pro-cessing , vol. 25, no. 3, pp. 1136–1151, March 2016.[6] T. Imbiriba, J. C. M. Bermudez, and C. Richard, “Band selectionfor nonlinear unmixing of hyperspectral images as a maximal cliqueproblem,”
IEEE Transactions on Image Processing , vol. 26, no. 5, pp.2179–2191, May 2017.[7] A. Zare and K. C. Ho, “Endmember variability in hyperspectral anal-ysis: Addressing spectral variability during spectral unmixing,”
SignalProcessing Magazine, IEEE , vol. 31, p. 95, January 2014.[8] L. Drumetz, J. Chanussot, and C. Jutten, “Variability of the endmembersin spectral unmixing: recent advances,” in ,Los Angeles, USA, 2016.[9] R. A. Borsoi, T. Imbiriba, and J. C. M. Bermudez, “Super-resolutionfor hyperspectral and multispectral image fusion accounting for seasonalspectral variability,”
IEEE Transactions on Image Processing , vol. 29,no. 1, pp. 116–127, 2020.[10] P.-A. Thouvenin, N. Dobigeon, and J.-Y. Tourneret, “Hyperspectralunmixing with spectral variability using a perturbed linear mixingmodel,”
IEEE Trans. Signal Processing , vol. 64, no. 2, pp. 525–538,Feb. 2016.[11] L. Drumetz, M.-A. Veganzones, S. Henrot, R. Phlypo, J. Chanussot,and C. Jutten, “Blind hyperspectral unmixing using an extended linearmixing model to address spectral variability,”
IEEE Transactions onImage Processing , vol. 25, no. 8, pp. 3890–3905, 2016.[12] T. Imbiriba, R. A. Borsoi, and J. C. M. Bermudez, “Generalized linearmixing model accounting for endmember variability,” in
ICASSP, IEEEInternational Conference on Acoustics, Speech and Signal Processing ,April 2018, pp. 1862–1866.[13] A. Halimi, J. Bioucas-Dias, N. Dobigeon, G. S. Buller, and S. McLaugh-lin, “Fast hyperspectral unmixing in presence of nonlinearity or mismod-elling effects,”
IEEE Trans. Computational Imaging , vol. 3, no. 2, pp.146–159, April 2017.[14] D. Hong, N. Yokoya, J. Chanussot, and X. X. Zhu, “An augmented linearmixing model to address spectral variability for hyperspectral unmixing,”
IEEE Transactions on Image Processing , vol. 28, no. 4, pp. 1923–1938,2019.[15] M. Tao and X. Yuan, “Recovering low-rank and sparse components ofmatrices from incomplete and noisy observations,”
SIAM Journal onOptimization , vol. 21, no. 1, pp. 57–81, 2011.[16] M. Bousse, O. Debals, and L. De Lathauwer, “A tensor-based methodfor large-scale blind source separation using segmentation,”
IEEE Trans-actions on Signal Processing , vol. 65, no. 2, pp. 346–358, 2017.[17] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalex-akis, and C. Faloutsos, “Tensor decomposition for signal processing andmachine learning,”
IEEE Transactions on Signal Processing , vol. 65,no. 13, pp. 3551–3582, 2017.[18] M. K.-P. Ng, Q. Yuan, L. Yan, and J. Sun, “An adaptive weighted tensorcompletion method for the recovery of remote sensing images withmissing data,”
IEEE Transactions on Geoscience and Remote Sensing ,vol. 55, no. 6, pp. 3367–3381, 2017.[19] X. Zhang, G. Wen, and W. Dai, “A tensor decomposition-based anomalydetection algorithm for hyperspectral image,”
IEEE Transactions onGeoscience and Remote Sensing , vol. 54, no. 10, pp. 5801–5820, 2016.[20] X. Guo, X. Huang, L. Zhang, L. Zhang, A. Plaza, and J. A. Benedik-tsson, “Support tensor machines for classification of hyperspectralremote sensing imagery,”
IEEE Transactions on Geoscience and RemoteSensing , vol. 54, no. 6, pp. 3248–3264, 2016.[21] S. Yang, M. Wang, P. Li, L. Jin, B. Wu, and L. Jiao, “Compressivehyperspectral imaging via sparse tensor and nonlinear compressed sens-ing,”
IEEE Transactions on Geoscience and Remote Sensing , vol. 53,no. 11, pp. 5943–5957, 2015.[22] L. Zhang, L. Zhang, D. Tao, and X. Huang, “Tensor discrimina-tive locality alignment for hyperspectral image spectral–spatial feature extraction,”
IEEE Transactions on Geoscience and Remote Sensing ,vol. 51, no. 1, pp. 242–256, 2013.[23] M. A. Veganzones, J. E. Cohen, R. C. Farias, J. Chanussot, andP. Comon, “Nonnegative tensor cp decomposition of hyperspectral data,”
IEEE Transactions on Geoscience and Remote Sensing , vol. 54, no. 5,pp. 2577–2588, 2016.[24] Y. Qian, F. Xiong, S. Zeng, J. Zhou, and Y. Y. Tang, “Matrix-vectornonnegative tensor factorization for blind unmixing of hyperspectral im-agery,”
IEEE Transactions on Geoscience and Remote Sensing , vol. 55,no. 3, pp. 1776–1792, 2017.[25] T. Imbiriba, R. A. Borsoi, and J. C. M. Bermudez, “A low-rank tensorregularization strategy for hyperspectral unmixing,”
IEEE StatisticalSignal Processing Workshop , pp. 373 – 377, 2018.[26] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinearsingular value decomposition,”
SIAM journal on Matrix Analysis andApplications , vol. 21, no. 4, pp. 1253–1278, 2000.[27] R. A. Borsoi, T. Imbiriba, Bermudez, and J. C. M., “A data dependentmultiscale model for hyperspectral unmixing with spectral variability,”
IEEE Transactions on Image Processing (Submitted) , 2020.[28] T. Uezato, M. Fauvel, and N. Dobigeon, “Hyperspectral unmixing withspectral variability using adaptive bundles and double sparsity,”
IEEETransactions on Geoscience and Remote Sensing , 2019.[29] L. Drumetz, T. R. Meyer, J. Chanussot, A. L. Bertozzi, and C. Jutten,“Hyperspectral image unmixing with endmember bundles and groupsparsity inducing mixed norms,”
IEEE Transactions on Image Process-ing , 2019.[30] R. A. Borsoi, T. Imbiriba, and J. C. M. Bermudez, “Deep generativeendmember modeling: An application to unsupervised spectral unmix-ing,”
IEEE Transactions on Computational Imaging , 2019, doi: 10.1109/TCI.2019.2948726.[31] A. Cichocki, D. Mandic, L. De Lathauwer, G. Zhou, Q. Zhao, C. Caiafa,and H. A. Phan, “Tensor decompositions for signal processing applica-tions: From two-way to multiway component analysis,”
IEEE SignalProcessing Magazine , vol. 32, no. 2, pp. 145–163, 2015.[32] J. H˚astad, “Tensor rank is np-complete,”
Journal of Algorithms , vol. 11,no. 4, pp. 644–654, 1990.[33] S. Gandy, B. Recht, and I. Yamada, “Tensor completion and low-n-ranktensor recovery via convex optimization,”
Inverse Problems , vol. 27,no. 2, p. 025010, 2011.[34] H. Fan, Y. Chen, Y. Guo, H. Zhang, and G. Kuang, “Hyperspectral imagerestoration using low-rank tensor recovery,”
IEEE Journal of SelectedTopics in Applied Earth Observations and Remote Sensing , vol. 10,no. 10, pp. 4589–4604, 2017.[35] S. Mei, J. Hou, J. Chen, L. P. Chau, and Q. Du, “Simultaneousspatial and spectral low-rank representation of hyperspectral images forclassification,”
IEEE Transactions on Geoscience and Remote Sensing ,vol. PP, no. 99, pp. 1–15, 2018.[36] L. Sorber, M. Van Barel, and L. De Lathauwer, “Optimization-based al-gorithms for tensor decompositions: Canonical polyadic decomposition,decomposition in rank-(Lr,Lr,1) terms, and a new generalization,”
SIAMJournal on Optimization , vol. 23, no. 2, pp. 695–720, 2013.[37] Q. Wang, Z. Qin, F. Nie, and X. Li, “Spectral embedded adaptive neigh-bors clustering,”
IEEE transactions on neural networks and learningsystems , no. 99, pp. 1–7, 2018.[38] Q. Wang, F. Zhang, and X. Li, “Optimal clustering framework forhyperspectral band selection,”
IEEE Transactions on Geoscience andRemote Sensing , no. 99, pp. 1–13, 2018.[39] T. Yokota, N. Lee, and A. Cichocki, “Robust multilinear tensor rankestimation using higher order singular value decomposition and infor-mation criteria,”
IEEE Transactions on Signal Processing , vol. 65, no. 5,pp. 1196–1206.[40] J. M. P. Nascimento and J. M. Bioucas-Dias, “Vertex Component Anal-ysis: A fast algorithm to unmix hyperspectral data,”
IEEE Transactionson Geoscience and Remote Sensing , vol. 43, no. 4, pp. 898–910, April2005.[41] R. N. Clark, G. A. Swayze, K. E. Livo, R. F. Kokaly, S. J. Sutley,J. B. Dalton, R. R. McDougal, and C. A. Gent, “Imaging spectroscopy:Earth and planetary remote sensing with the usgs tetracorder and expertsystems,”
Journal of Geophysical Research: Planets , vol. 108, no. E12,2003.[42] B. Hapke, “Bidirectional reflectance spectroscopy, 1, Theory,”
Journalof Geophysical Research , vol. 86, no. B4, pp. 3039–3054, 1981.[43] R. A. Borsoi, T. Imbiriba, and J. C. Moreira Bermudez, “Improvedhyperspectral unmixing with endmember variability parametrized usingan interpolated scaling tensor,” in , May 2019, pp.2177–2181. [44] R. A. Borsoi, T. Imbiriba, J. C. M. Bermudez, and C. Richard, “Afast multiscale spatial regularization for sparse hyperspectral unmixing,”
IEEE Geoscience and Remote Sensing Letters , vol. 16, no. 4, pp. 598–602, 2019.
Tales Imbiriba (S’14, M’17) received his Doctoratedegree from the Department of Electrical Engi-neering (DEE) of the Federal University of SantaCatarina (UFSC), Florian´opolis, Brazil, in 2016. Heserved as a Postdoctoral Researcher at the DEE–UFSC and is currently a Postdoctoral Researcherat the ECE dept. of the Northeastern University,Boston, MA, USA. His research interests include au-dio and image processing, pattern recognition, kernelmethods, adaptive filtering, and Bayesian Inference.
Ricardo Augusto Borsoi (S’18) received the MScdegree in electrical engineering from Federal Uni-versity of Santa Catarina (UFSC), Florian´opolis,Brazil, in 2016. He is currently working towards hisdoctoral degree at Universit´e Cˆote d’Azur (OCA)and at UFSC. His research interests include imageprocessing, tensor decomposition, and hyperspectralimage analysis.