Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization

0.71 MB·

by  Cédric Févotte

1
Nonlinear hyperspectral unmixing with
robust nonnegative matrix factorization
Ce´dric Fe´votte and Nicolas Dobigeon
Abstract—This paper introduces a robust mixing model to using bilinear models [8], [9]. As explained in [10], several
describehyperspectraldataresultingfromthemixtureofseveral bilinear models have been proposed [11]–[13], and they
pure spectral signatures. This new model not only general-
mainly differ by the constraints imposed on the nonlinearity
izes the commonly used linear mixing model, but also allows
4 for possible nonlinear effects to be easily handled, relying on term. The linear-quadratic model advocated in [14] also
1 mild assumptions regarding these nonlinearities. The standard incorporates pairwise interactions between the endmembers
0 nonnegativity and sum-to-one constraints inherent to spectral components. Conversely, to approximate a large range of
2 unmixing are coupled with a group-sparse constraint imposed second-order nonlinearities, Altmann et al. [15] introduce a
r on the nonlinearity component. This results in a new form polynomial post-nonlinear model that has demonstrated its
a of robust nonnegative matrix factorization. The data fidelity
M term is expressed as a β-divergence, a continuous family of ability to describe most of the nonlinear effects, in particular
dissimilaritymeasuresthattakesthesquaredEuclideandistance in vegetated areas [16]. A common feature of these models

6 andthegeneralizedKullback-Leiblerdivergenceasspecialcases. is that they all consist in including a supplementary additive
The penalized objective is minimized with a block-coordinate term to the standard LMM, accounting for the nonlinearities.

] descent that involves majorization-minimization updates. Sim- One major drawback of these models, however, is that they
E ulation results obtained on synthetic and real data show that
require to choose a specific form of nonlinearity, and this can
M the proposed strategy competes with state-of-the-art linear and
nonlinear unmixing methods. be limiting in practice.
.
t
a Index Terms—Hyperspectral imagery, nonlinear unmixing,
In this paper, a new so-called robust LMM (rLMM) is
t robust nonnegative matrix factorization, group-sparsity.
s proposed. Similarly to the nonlinear models detailed above, it
[
is built on the standard LMM and includes a supplementary

2 I. INTRODUCTION additive term that accounts for nonlinear effects. However,
v Spectral unmixing (SU) is an issue of prime interest it does not require to specify an analytical form of the
9
when analyzing hyperspectral data since it provides a nonlinearity. Instead, nonlinearities are merely treated as
4
6 comprehensive and meaningful description of the collected outliers. The primary motivation is that the LMM can be
5 measurements in various application fields including remote considered as a valid model to describe the majority of the
1. sensing [1], planetology [2], food monitoring [3] or spectro- pixels in a remotely sensed image and, as a consequence,
0 microscopy [4]. SU consists in decomposing P multi-band only a sparse number of pixels are affected by nonlinearities.
4 observations Y = [y ,...,y ] into a collection of K As such, one of the contributions reported in this article
1 P
1 individual spectra M = [m ,...,m ], called endmembers, consists in decomposing the L×P matrix of the multi-band
1 K
v: and estimating their relative proportions (or abundances) observations as Y ≈ MA + R, where R is a sparse
i A = [a ,...,a ] in each observation [5], [6]. Most of the (and nonnegative) residual term accounting for outliers
X 1 P
hyperspectral unmixing algorithms proposed in the signal (i.e., nonlinear effects). To reflect the assumption that the
ar & image processing and geoscience literatures rely on the LMM holds for most of the observed pixels, the sparsity
commonly admitted linear mixing model (LMM), Y ≈MA. constraint is imposed at the group-level, i.e., a column of
Indeed, LMM provides a good approximation of the physical R will be assumed to be either entirely zero or not. The
process underlying the observations and has resulted in proposed decomposition relates to robust nonnegative matrix
interesting results for most applications. However, for several factorization (rNMF) as will be explained in more details in
specific applications, LMM may be inaccurate and other the sequel of the paper.
nonlinear models need to be advocated [7]. For instance, in
remotely sensed images composed of vegetation (e.g., trees), The article is organized as follows. The rLMM is intro-
interactions of photons with multiple components of the duced in more details in Section II. Section III describes
scene lead to nonlinear effects that can be taken into account a block-coordinate descent algorithm for rLMM estimation.
Experimental results obtained on synthetic data are reported
ThisworkissupportedbytheESTOMATPEPSProjectsupportedbyCNRS in Section IV. Two real hyperspectral images are investigated
andbytheHypanemaANRProjectn◦ANR-12-BS03-003.
in Section V. Section VI concludes. This article extends our
C.Fe´votteiswithLaboratoireLagrange(CNRS,Universite´deNiceSophia
Antipolis and Observatoire de la Coˆte d’Azur), Parc Valrose, 06108 Nice preliminary conference paper [17] in a significant way. We
cedex2,France.(e-mail:[email protected]). here generalize the use of the squared Euclidean distance
N. Dobigeon is with University of Toulouse, IRIT/INP-ENSEEIHT, 2
considered in [17] to the more general β-divergence. Addi-
rue Camichel, BP 7122, 31071 Toulouse cedex 7, France. (e-mail: nico-
[email protected]). tionally, we show how some of the multiplicative updates
2
obtained heuristically in [17] can be rigorously obtained via Eq. (5) defines a robust NMF problem. Robust NMF is a
majorization-minimization. We also describe a rule of thumb nonnegativevariantofrobustPCA[18]whichhasappearedin
forchoosingthevalueofthepenaltyweightefficiently.Finally, different forms in the literature. In [19], the outlier term R is
we provide extended experimental results on synthetical and nonnegativeandpenalizedbythe(cid:96) norm.In[20]and[21],R
1
real data. isreal-valuedandpenalizedby(cid:96) and(cid:96) norms,respectively.
1 1,2
In [22], the (cid:96) norm of (Y−MA) is minimized (noise free
2,1
II. ROBUSTLINEARMIXINGMODEL scenario). A so-called robust nonnegative matrix factorization
A. Model design approach was introduced for the reconstruction of reflectance
spectra in [23]; however the term “robust” there refers to a
The proposed rLMM is described by
different feature, namely the use of a data-fitting term (the
(cid:88)K hypersurface cost function) that is less sensitive to outlier
y ≈ a m +r , (1)
p kp k p observations than the traditional squared Euclidean distance,
k=1 forthecomputationofaregularNMFY ≈MA.Notefinally
where y = [y ,...,y ]T denotes the pth pixel spectrum that other articles that have addressed hyperspectral unmixing
p 1p Lp
observed in L spectral bands, m = [m ,...,m ]T de- with regular NMF (i.e., in the standard linear model), e.g.,
k 1k Lk
notes the kth endmember spectrum, a = [a ,...,a ]T [24]–[27].
p 1p Kp
denotes the abundances associated with the pth pixel and To the best of our knowledge, the formulation of robust
r =[r ,...,r ]T denotes the outlier term (accounting for NMF described by Eq. (5), where R is nonnegative and
p 1p Lp
nonlinearities). The matrix formulation of Eq. (1) is given by penalized by the (cid:96)2,1 norm (and where the abundances sum
to one), is entirely novel. Furthermore, and most importantly,
Y ≈MA+R. (2)
previouswork[19]–[22]hasonlyconsideredrobustNMFwith
The approximation symbol in Eqs. (1) and (2) underlies the the squared Euclidean distance, i.e., D(Y|MA+R)=(cid:107)Y−
minimization of a measure of dissimilarity D(Y|MA+R), MA−R(cid:107)2 while we give here a more general formulation
2
the β-divergence, that will be introduced in Section II-C. based on the β-divergence, defined next.
The matrices Y, M and A are nonnegative by nature and
we take the abundance coefficients to sum to one, i.e., C. The β-divergence
(cid:40) K (cid:41) We take the measure of fit in Eq. (5) to be such that
ap ∈SK d=ef a∈RK(cid:12)(cid:12) ak ≥0, (cid:88)ak =1 , (3) D(A|B)=(cid:88)d(a |b ), (7)
ij ij
k=1
ij
as commonly assumed in most hyperspectral data models.
where d(x|y) is the β-divergence between positive scalars x
In this work, we assume the nonlinear component r to be
p and y. The β-divergence was introduced in various forms
nonnegative as well, like in the bilinear models of [11]–[13]
in [28]–[30] and has become a standard measure of fit in
and the polynomial model with constructive interferences of
NMF, see, e.g., [31]–[34]. In this paper we use the following
[15].Asdiscussedintheintroduction,weexpectr tobeoften
p definition:
zero, i.e., pixels to follow the standard LMM in general. For

pixels where the LMM assumption fails, nonlinearities will  β(xββ−1) + yββ − xβy−β−11 β ∈R\{0,1}
def
become “active” and r will become nonzero. This amounts d (x|y) = x logx −x+y β =1 . (8)
p β y
to say that the energy vector  x −logx −1 β =0
y y
e=[(cid:107)r (cid:107) ,...,(cid:107)r (cid:107) ] (4) The limiting cases β = 0 and β = 1 correspond to the
1 2 P 2
Itakura-Saito and Kullback-Leibler divergences, respectively.
issparse.InEq.(4),(cid:107)·(cid:107) denotestheEuclideannormdefined
by (cid:107)x(cid:107) = (cid:112)(cid:80) x2. S2parsity can routinely be enforced by The squared Euclidean distance is obtained for β = 2, i.e.,
2 k k d (x|y)=(x−y)2/2.Theparameterβ essentiallygoverns
(cid:96) -regularization, as done next. β=2
1 the assumed statistics of the observation noise and can either
befixedorlearntfromtrainingdatabycross-validation.Under
B. Objective function
certain assumptions, the β-divergence can be mapped to a
In light of previous section, our objective is to solve the log-likelihoodfunctionfortheTweediedistribution[28],[35],
minimization problem defined by [36], parametrized with respect to its mean. In our setting,
this translates into E[Y|M,A,R]=MA+R. In particular,
min J(M,A,R)=D(Y|MA+R)+λ(cid:107)R(cid:107)
2,1 the values β = 0,1,2 underlie multiplicative Gamma noise,
M,A,R
PoissonnoiseandGaussianadditivenoise,respectively.Theβ-
s.t. M≥0, A≥0, R≥0 and (cid:107)a (cid:107) =1, (5)
p 1 divergence offers a continuum of noise statistics that interpo-
where λ is a nonnegative penalty weight, A ≥ 0 denotes latesbetweenthesethreespecificcases.Anoteworthyproperty
(cid:80)
nonnegativity of the coefficients of A, (cid:107)x(cid:107)1 = k|xk| and of the β-divergence is its behavior with respect to scale.
(cid:107)·(cid:107)2,1 is the so-called (cid:96)2,1-norm defined by Indeed, let λ ∈ R+, then we have dβ(λx|λy) = λβdβ(x|y).
Thismeansthatthedata-fittingtermwillrelymoreheavilyon
P
(cid:88)
(cid:107)R(cid:107) =(cid:107)e(cid:107) = (cid:107)r (cid:107) . (6) large (resp., small) coefficients in Y for β >0 (resp., β <0);
2,1 1 p 2
see a more detailed discussion in [37].
p=1
3
TABLEI
III. BLOCK-COORDINATEDESCENTALGORITHM DIFFERENTIABLECONVEX-CONCAVEDECOMPOSITIONOFTHE
InordertosolvetherNMFminimizationproblemdefinedat β-DIVERGENCEANDMMUPDATEEXPONENTS.
Eq.(5),wepresentaniterativeblock-coordinatedescentalgo-
(cid:94) (cid:95)
d(x|y) d(x|y) γ(β) ξ(β)
rithmthatupdateseachoftheparametersM,AandRinturn.
Eachparameterisupdatedconditionallyuponthecurrentvalue β<1andβ(cid:54)=0 − 1 xyβ−1 1yβ 1 1
β−1 β 2−β 3−β
of the other parameters and such that the objective function is β=0 xy−1 logy 1 1
decreased. This is the updating scheme employed by virtually 2 3
1≤β≤2 d(x|y) 0 1 1
all NMF algorithms. Unfortunately, given the non-convexity 3−β
of the objective function J(M,A,R), this strategy can return β>2 β1yβ −β−11xyβ−1 β−11 β−11
local solutions and proper initialization is required. This will
be addressed in Section IV. The updates of the parameters
are described next. In short, the parameters M and R are that A and R are here treated as constants). Then, define
updated via majorization-minimization (MM). The parameter for k = 1,...,K, λ˜ = m˜ a /y˜ and for k = K +1,
lkp lk kp lp
A is updated using a heuristic scheme that has proven to λ˜ = r /y˜ . By construction, we have (cid:80)K+1λ˜ = 1.
lkp lp lp k=1 lkp
work well in the literature. All the updates turn out to be Then, using definition of convexity, we have
“multiplicative”, i.e., such that the new update is obtained by
(cid:32) (cid:33)
term-to-termmultiplyingthepreviousupdatebyanonnegative C(cid:94)(M)=(cid:88)(cid:94)d y |(cid:88)λ˜ mlkakp +λ˜ rlp
matrix,henceautomaticallypreservingthenonnegativityofthe lp lkp λ˜ l(K+1)pλ˜
estimatesthroughiterations.Theresultingalgorithmhaslinear lp k lkp l(K+1)p
complexity O(LKP) (in flops) per iteration. ≤(cid:88)(cid:34)(cid:88)K λ˜ (cid:94)d(cid:32)y |mlkakp(cid:33)
lkp lp λ˜
lp k=1 lkp
A. Update of the endmember spectra M (cid:32) (cid:33)(cid:35)
Updating M given the current values of A and R involves + λ˜l(K+1)p(cid:94)d ylp|λ˜ rlp
solving the following minimization problem l(K+1)p
(cid:34) K (cid:18) (cid:19) (cid:35)
min C(M)=D(Y|MA+R) s.t. M≥0. (9) =(cid:88) (cid:88)m˜lkakp(cid:94)d y |y˜ mlk + rlp(cid:94)d(y |y˜ )
M y˜lp lp lpm˜lk y˜lp lp lp
lp k=1
When R = 0, this problem boils down to updating the
d=efG(cid:94)(M|M˜) (12)
dictionary matrix in NMF with the β-divergence. MM al-
gorithms have been designed for that purpose in [33], [34].
The auxiliary function essentially “breaks” the sum over k in
In this section, we extend the MM approach to the case
where R ≥ 0. Denote by M˜ the estimate of M at current Eq. (10) to make the optimization over M separable w.r.t its
entries m .
iteration. The first step of MM consists in building an upper lk
(cid:95)
bound G(M|M˜) of C(M) which is tight for M = M˜, i.e., Thanks to its concavity, C(M) can be majorized by a first-
C(M) ≤ G(M|M˜) for all M and C(M˜) = G(M˜|M˜). The order approximation at M = M˜ (the tangent inequality),
second step consists in minimizing the bound with respect leading to
to (w.r.t) M, producing a valid descent algorithm. Indeed, at
iteration i + 1, it holds by construction that C(M(i+1)) ≤ G(cid:95)(M|M˜)=C(cid:95)(M˜)+(cid:88)(cid:95)d(cid:48)(y |y˜ )(cid:88)a (m −m˜ ),
lp lp kp lk lk
G(M(i+1)|M(i)) ≤ G(M(i)|M(i)) = C(M(i)). The bound lp k
G(M|M˜) will be referred to as auxiliary function. (13)
The auxiliary function obtained in [34] relies on a convex-
(cid:95) (cid:95)
concave decomposition of d(x|y) and we follow a similar where d(cid:48)(x|y) denotes the derivative of d(x|y) w.r.t y.
approach here. The β-divergence can always be decomposed AnupperboundtoC(M)isfinallyobtained(uptoconstant
as d(x|y)=(cid:94)d(x|y)+(cid:95)d(x|y)+cst where (cid:94)d(x|y) and (cid:95)d(x|y) terms) by adding G(cid:94)(M|M˜) and G(cid:95)(M|M˜). Skipping details
are respectively convex and concave functions of y and cst is forbrevity,theresulting functioncanbeminimizedinclosed-
constant w.r.t y. Such a decomposition is not unique; we give form w.r.t M˜, resulting in the following update
a “natural” decomposition in Table I. It follows that C(M)
can be decomposed as the sum of a convex term C(cid:94)(M), a (cid:32)(cid:80) a y y˜β−2(cid:33)γ(β)
concave term C(cid:95)(M) and a constant term, such that mlk =m˜lk (cid:80)p kap ly˜pβ−lp1 , (14)
p kp lp
(cid:94) (cid:88)(cid:94) (cid:88)
C(M)= d(y | m a +r ), (10)
lp lk kp lp
where γ(β) is given in Table I and we recall that y˜ =
lp k lp
(cid:80)
C(cid:95)(M)=(cid:88)(cid:95)d(y |(cid:88)m a +r ). (11) km˜lkakp+rlp isthedataapproximationatcurrentiteration.
lp lk kp lp Note that we observed in practice that setting γ(β) = 1
lp k
for every value of β still reduced the objective function at
(cid:94)
Fromthere,C(M)canbemajorizedusingJensen’sinequal- every iteration and produced faster convergence. This may be
(cid:80)
ity, as follows. Let us denote y˜ = m˜ a +r the data interpreted as over-relaxation of the MM update, see [34] for
lp k lk kp lp
approximation formed with the current iterate M˜ (and recall further discussion on this subject.
4
B. Update of the outlier term R C. Update of the abundances A
Updating R given the current values of M and A involves Updating A given the current values of M and R involves
solving the following minimization problem solving the following minimization problem
minC(A)=D(Y|MA+R) s.t. A≥0and∀p,(cid:107)a (cid:107) =1.
p 1
minC(R)=D(Y|MA+R)+λ(cid:107)R(cid:107)2,1 s.t. R≥0. (15) A
R (19)
The data-fitting term may be majorized using a convex- Thesum-to-oneconstraintontheabundancesinducesanextra
concave decomposition of D(Y|MA + R) exactly as we difficulty as compared to the optimization problems involved
did in Section III-A. Denote by R˜ the current update of R, by the updates of M and R. In some cases such a constraint
(cid:80)
s = m a = [MA] the low-rank component and can be handled using Lagrange multipliers, but this approach
lp k lk kp lp
y˜ =s +r˜ thecurrentdataapproximation.1Then,applying does not succeed in our setting, except in the special case
lp lp lp
the Jensen and tangent inequalities to the convex and concave β = 1, corresponding to the generalized Kullback-Leibler
parts, respectively, we obtain divergence. We hence resort to another common approach
based on a change of variable. We introduce the variable U
(cid:20) (cid:21)
D(Y|MA+R)≤(cid:88) r˜lp(cid:94)d(y |y˜ rlp)+ slp(cid:94)d(y |y˜ ) to be a nonnegative matrix of dimension K×P and set
lp y˜lp lp lpr˜lp y˜lp lp lp a = ukp . (20)
(cid:88) (cid:95) (cid:95) kp (cid:107)u (cid:107)
+ [d(y |y˜ )+ d(cid:48)(y |y˜ )(r −r˜ )]. p 1
lp lp lp lp lp lp
The optimization problem of Eq. (19) is turned into the new
lp
(16) optimization problem
(cid:18) (cid:20) (cid:21) (cid:19)
u u
Denote by F(R|R˜) the right-hand side of Eq. (16). An minC(U)=D Y |M 1 ,..., P +R
auxiliary function for C(R) may simply be obtained as U (cid:107)u1(cid:107)1 (cid:107)uP(cid:107)1
G(R|R˜) = F(R|R˜) + λ(cid:107)R(cid:107) . However, this specific s.t. U≥0 (21)
2,1
auxiliary function is not amenable to optimization w.r.t R (no which is free from the sum-to-one constraint. This approach
closed-form solution). Hence, the first step of our strategy is has been used for NMF in [38]. Unfortunately, we were not
tomajorizethepenaltyfunction(cid:107)R(cid:107)2,1 aswell.Byconcavity able to produce an auxiliary function for the new objective
of the square-root function, we may write function in (21) – in particular because it can no longer be
easilydecomposedasaconvexpartandconcavepart.Instead,
(cid:107)R(cid:107) ≤ 1(cid:88)(cid:18)(cid:107)rp(cid:107)22 +(cid:107)˜r (cid:107) (cid:19). (17) weresorttoaheuristiccommonlyusedinNMF,see,e.g.,[37],
2,1 2 (cid:107)˜r (cid:107) p 2
p 2 [39], as follows. As it appears, the gradient of C(U) can be
p
expressed as the difference of two nonnegative functions such
(cid:113)
Equation (17) essentially replaces (cid:80) r2 by a quadratic that
l lp
tight upper-bound that involves (cid:80) r2, with the effect of ∇ C(U)=∇+ C(U)−∇− C(U). (22)
l lp ukp ukp ukp
decouplingthespectralbandsfromwithinthesquareroot.Un-
The heuristic algorithm simply writes
fortunately,theresultingauxiliaryfunctionisnotyetamenable
to optimization. As such, from here our approach closely ∇− C(U˜)
ofoflldoewgsre[e36β].,Fmoratβch>ing2,twhee mmaoynommaijaolrizoefrhl2pigbhyesat mdeognroemeiianl ukp =u˜kp∇u+ukkppC(U˜). (23)
F(R|R˜) (see Table I). For β ≤ 2, the reverse is done: the It ensures nonnegativity of the parameter updates provided
leading monomial in F(R|R˜) is now of degree lower than 2 initializationwithanonnegativevalue,andproducesadescent
algorithm in the sense that u is updated towards left (resp.,
and may as such be majorized by a quadratic term, matching kp
right) when the gradient is positive (resp., negative). The
thequadraticupperboundofthepenaltyfunction;seeSection
algorithm was found experimentally to decrease the value of
4.1in[36]formoredetails.Thisstrategyleadstothefollowing
the objective function at each iteration for every of the many
update
(cid:80)
values of β that we tried. Denoting s˜ = m a˜ and
lp k lk kp
 ξ(β) y˜ =s˜ +r , the update is found to be
y y˜β−2 lp lp lp
rlp =r˜lpy˜β−1lp+lλp r˜lp  , (18) u =u˜ (cid:80)l(mlkylpy˜lβp−2+s˜lpy˜lβp−1). (24)
lp (cid:107)˜rp(cid:107)2 kp kp(cid:80) (m y˜β−1+s˜ y y˜β−2)
l lk lp lp lp lp
where ξ(β) is the exponent given in Table I. Again, we
The update for A is then simply a =u /(cid:107)u (cid:107) .
kp kp p 1
observed in practice that setting ξ(β) = 1 for every value As it turns out, the updates (14), (18) and (23) can be
of β still reduced the objective function at every iteration and implemented in matrix form, as shown in Algorithm 1, which
produced faster convergence. recapitulates the overall procedure. In Algorithm 1, all op-
erators preceded by a dot ‘·’ are entrywise MATLAB-like
Sec1TtiohneIsIaI-mAeanndotfaotiroyn˜lpy˜l=p(cid:80)iskumseldkafkopr+y˜lr˜plp=inS(cid:80)ectkiom˜nlIkIIa-kBp.O+urrilnpteinnt owpeelrl.atAiodndsitaionndalflrya,ct1ion badrsensohtaelsl tbheetMake×n tNermm-taot-rtiexrmwiaths
is too avoid the use of too many notations and the definition of y˜lp should M,N
beclearfromcontext(i.e.,whichparameterupdateisconsidered). coefficients equal to 1.
5
Algorithm 1 Group robust NMF Let us now assume an unspecified independent prior model
Initialize M, A and R for MA but such that E[[MA] ] = ρ. Denoting by µˆ =
lp
Set convergence tolerance parameter ‘tol’ (LP)−1(cid:80)y the empirical data expectation, our approach
lp
Set value of λ consists in matching µˆ with E[[MA] ]+E[r ], leading to
lp lp
S=MA
C
Yˆ =S+R λˆ = . (26)
µˆ−ρ
while err ≥ tol do
% Update outlier term R We insist that the latter expression only provides a handy
gross estimate of λ that comes with no statistical guarantee.
(cid:34) (cid:35) In particular the estimate of λ is extremely dependent on
Y.Yˆ.(β−2)
R←R. ρ, the prior expectation of [MA]lp. However, because ρ is
Yˆ.(β−1)+λRdiag[(cid:107)r1(cid:107)1,...,(cid:107)rP(cid:107)1]−1 lower bounded by 0, the estimate of λ is lower bounded by
Yˆ ←S+R λ0 = C/µˆ, corresponding to a plausible minimum degree of
sparsity. We used λ = λ in the evaluations below and this
0
was found to provide satisfactory results.
% Update abundances A
MT(Y.Yˆ.(β−2))+1 (S.Yˆ.(β−1)) IV. EXPERIMENTSWITHSYNTHETICDATA
A←A. K,L
MT(Yˆ.(β−1))+1 (S.Y.Yˆ.(β−2)) In this section we evaluate the relevance of the rLMM
K,L
proposed in Section II and the accuracy of the corresponding
S←MA
rNMF algorithm described in Section III using synthetic data.
Yˆ ←S+R
A. Data generation
% Update endmembers M Four 64×64-pixel images composed of K = 3 or 6 pure
(cid:34)(Y.Yˆ.(β−2))AT(cid:35) spectral components have been generated according to four
M←M. differentlinearandnonlinearmodels.Theendmemberspectra
(Yˆ.(β−1))AT
have been extracted from the spectral library provided with
S←MA theENVIsoftware[41].Thefirstimage,denotedasI ,is
LMM
Yˆ ←S+R composedofpixelsfollowingthestandardLMM(nononlinear
component)
K
(cid:88)
Computetheobjectivefunctionrelativedecrease‘err’(or yp = akpmk+np, (27)
any other convergence criterion). k=1
end while with a ∈ SK. The three other images, denoted I , I
p NM FM
and I , are generated as follows. Three fourths of the
GBM
image pixels are generated according to the LMM in (27) and
D. Setting the value of λ the remaining fourth is generated according to a model that
featuresnonlinearcomponentinteractions.Moreprecisely,the
The hyperparameter λ controls the trade-off between the
latter pixels are generated according to:
data-fittingtermD(Y|MA+R)andthepenaltyterm(cid:107)R(cid:107) .
2,1
Setting the “right” value of λ is a difficult task, like in any • the Nascimento model (NM) [11]
otherso-calledvariationalapproachthatinvolvesaregulariza- K K−1 K
(cid:88) (cid:88) (cid:88)
tion term. We describe in this paragraph a rule of thumb for y = a m + b m (cid:12)m +n ,
p kp k ip i j p
choosing λ in a plausible range of values. Our approach is k=1 i=1 j=i+1
based on the method of moments. It consists in interpreting
with
theobjectivefunction(5)asajointlikelihoodandinmatching (cid:20) (cid:21)
a
the empirical mean of the data with its prior expectation in p ∈S2K−1 (28)
b
the statistical model. As mentioned in Section II-C, the β- p
divergence underlies a Tweedie data distribution such that and b =[b ,...,b ]T,
p 1p (K−1)p
E[Y|MA+R]=MA+R. The term λ(cid:107)R(cid:107)2,1 can be seen • the Fan bilinear model (FM) [12]
a log-prior term. Using some results from [40], the corre-
K K−1 K
spondingpriordistributionp(r )foreachcolumnofRcanbe (cid:88) (cid:88) (cid:88)
p y = a m + a a m (cid:12)m +n ,
p kp k ip jp i j p
obtained as a scale mixture of conditionally independent half-
k=1 i=1 j=i+1
Normal distributions, with a Gamma distribution assigned to
the scale parameter. In particular, the expectation of rpl under with ap ∈SK,
this prior can be found to be • the generalized bilinear model (GBM) [15]
K K−1 K
2 Γ(K/2+1) 1 def C (cid:88) (cid:88) (cid:88)
E[rlp]= √πΓ(K/2+1/2)λ = λ. (25) yp = akpmk+ γijpaipajpmi(cid:12)mj+np,
k=1 i=1 j=i+1
6
with a ∈SK and where the nonlinear coefficient γ ∈ the relative difference between two successive values of the
p ijp
(0,1) adjust the bilinear interaction between the ith and objective function fell under 10−5.
jth endmembers in the pth pixel.
In the models introduced above, mi (cid:12)mj stands for the C. Performance measures
termwise (Hadamard) product.
Theperformanceoftheunmixingalgorithmswasevaluated
in terms of endmember estimation accuracy using the average
In a first experiment, the four images I , I , I
LMM NM FM spectral angle mapper (aSAM)
and I have been generated by drawing the abundance
GBM
acnodeffiucniiefonrtsmalypf(roorm[athTpeibrTpad]Tmiisnsibthlee sceatseSKof(oIrNMS2)Kr−an1)d.oImnlya aSAM(M)= K1 (cid:88)K acos(cid:18)(cid:107)(cid:104)mmk(cid:107),(cid:107)mˆmˆk(cid:105)(cid:107)(cid:19)
k k
second experiment, we wanted to evaluate the robustness of k=1
thealgorithmsw.r.ttheabsenceofpurepixelsintheimagesto and abundance estimation accuracy using the global mean
be unmixed. To do so, we imposed a cutoff to the abundance square error (GMSE)
coefficientsthatremovespurepixelsfromtheobservations.As
P
such, in this case the abundances have been uniformly drawn GMSE2(A)= 1 (cid:88)(cid:107)a −aˆ (cid:107)2.
over a truncated version of the set defined by (3), namely KP p p
p=1
(cid:40) K (cid:41)
SK0.9 = a∈RK(cid:12)(cid:12) ak ≥0, (cid:88)ak ≤0.9 . (29) D. Results and discussion
k=1 The performance measures returned by the unmixing meth-
Finally, in the two experiments the interaction coefficients ods are reported in Table II. First, the aSAM values show
γ appearing in the GBM have been uniformly drawn over that the proposed rNMF algorithm competes favorably with
ijp
the set (0,1) and the additive noise n was chosen white the two considered state-of-the-art endmember extraction al-
p
Gaussianwithsignal-to-noiseratioSNR=30dB,whichisan gorithms, namely VCA and Heylen’s algorithm. Initialized by
admissible value for most of the real imaging spectrometers. these algorithms, it almost always improves the endmember
estimation accuracy, with or without pure pixels. Similarly,
when analyzing the GMSE related to abundance estimation,
B. Compared methods
theseresultsdemonstratetheflexibilityoftherLMMtomodel
ThefourimageshavebeenunmixedusingrNMFandstate- observations coming from various scenarios. More generally,
of-the-art algorithms specially designed for the considered these results demonstrate the ability of the rLMM-based un-
models. The state-of-the-art algorithms are two-steps; the mixingtechniquetomitigateseveralkindsofnonlineareffects
endmember matrix M is estimated in a first step, and then whilepreservinggoodestimationperformancewhenanalyzing
the abundance matrix A is estimated in a second step, given only linear mixtures.
the endmember estimates (in a so-called “inversion” step). In
contrast, rNMF performs a joint estimation of M and A (and V. EXPERIMENTSWITHREALDATA
R).
InthissectionweapplyrNMFtorealhyperspectraldatasets
We considered vertex component analysis (VCA) [42] cou-
and discuss the results.
pled with fully constrained least squares (FCLS) [43]. VCA
and FCLS are standard endmember extraction and inversion
A. Description of the datasets
methods designed for the LMM. Besides, we considered the
nonlinear endmember extraction technique proposed in [44], We consider two real hyperspectral images that have been
denoted as Heylen’s algorithm in what follows, coupled with chosenbecauseofavailabilityofpartialgroundtruth.Thefirst
fourdifferentinversionmethodsdesignedforvariousnonlinear image was acquired over Moffett Field, CA, in 1997, by the
models, namely the NM, FM, GBM and the very flexible theAirborneVisibleInfraredImagingSpectrometer(AVIRIS)
polynomial post-nonlinear mixing model (PPNM) [15]. NM [46]. Water absorption bands have been removed from the
inversion is also achieved with FCLS since the NM can 224spectralbands,leadingtoL=189spectralbandsranging
be interpreted as a linear mixture of an extended set of from0.4µmto2.5µmwithanominalbandwidthof10nm.The
endmembers[11].FMinversionisachievedwiththealgorithm sceneofinterest,ofsizeof50×50pixels,consistsofapartof
detailed in [12], which exploits a first-order Taylor series lake and a coastal area composed of soil and vegetation. This
expansion of the nonlinear term. GBM inversion is achieved datasethasbeenpreviouslystudiedin[13],[47]and,thus,the
with the gradient descent algorithm from [45]. Finally, PPNM unmixingresultsobtainedinthecurrentworkcanbecompared
inversionisaddressedwiththesubgradient-basedoptimization tothosereportedintheselaterreferences.Thisdatasetwillbe
scheme from [15]. referred to as the “Moffett” image in the following.
rNMF is applied with β = 2 (reflecting the Gaussian The second considered dataset was acquired by the Hyspex
additive noise used in the data generation and for fair com- hyperspectral scanner over Villelongue, France, in 2010. The
parison with the other methods that rely on this assumption sensed spectral domain consists of L = 160 spectral bands
as well) and λ = λ . We considered initializations by either ranging from 0.4µm to 1.0µm. This image consists of a
0
VCA or Heylen’s algorithm. Convergence was stopped when forestedareawhere12vegetationspecieshavebeenidentified,
7
TABLEII
ESTIMATIONPERFORMANCEINTERMOFASAM(M)(×10−3)ANDGMSE2(A)(×10−3).BESTSCORESAPPEARINBLUEBOLDFACEANDSECOND
BESTSCORESAPPEARINBLUE.RNMFISINITIALIZEDBYEITHERVCAORHEYLEN’SMETHOD,ASSTATEDBETWEENBRACKETS.REFERTOTEXTFOR
OTHERDETAILS.
aSAM(M) GMSE2(A)
rNMF rNMF VCA Heylen rNMF rNMF
VCA Heylen
(VCA) (Heylen) +FCLS +NM +FM +GBM +PPNM (VCA) (Heylen)
ILMM 9.71 33.86 7.65 6.78 0.10 0.06 4.32 0.07 0.07 0.07 0.04
s =3 INM 12.74 170.50 202.05 256.71 20.78 73.79 78.34 96.23 91.85 54.68 82.80
xel R IFM 7.12 102.26 7.55 29.61 0.84 13.54 38.17 13.28 14.83 0.71 1.57
pi I 8.26 33.10 5.69 5.02 0.25 0.23 4.01 0.07 0.08 0.20 0.20
e GBM
ur
p ILMM 72.33 86.10 23.03 24.92 2.74 2.90 26.03 2.84 2.76 0.27 0.77
w/o =6 INM 175.97 249.60 156.81 238.68 53.40 49.00 75.31 74.87 70.38 32.12 32.46
R IFM 93.28 187.55 144.74 259.01 13.67 34.71 68.13 39.01 34.51 10.41 19.86
I 79.24 106.52 83.08 85.49 4.55 5.73 32.21 5.50 4.31 4.95 4.57
GBM
ILMM 46.45 64.89 12.37 11.91 1.89 2.05 11.85 2.07 2.06 0.16 0.14
s =3 INM 46.40 176.29 189.76 254.58 23.94 66.17 67.18 91.05 92.78 52.77 68.75
xel R IFM 52.77 214.15 9.26 239.32 4.03 107.59 104.75 96.43 115.65 0.79 20.07
pi I 48.18 58.58 9.49 9.31 2.66 2.48 5.85 1.77 1.58 0.30 0.29
e GBM
ur
p ILMM 66.55 79.19 15.72 17.67 3.95 2.34 24.35 2.07 2.11 0.63 0.39
with =6 INM 224.81 112.82 101.07 265.25 63.85 16.88 68.20 39.56 41.66 30.93 28.94
R IFM 98.16 171.10 145.41 198.64 11.21 28.55 56.15 30.30 28.64 9.28 10.20
I 75.21 136.20 79.34 156.75 5.36 16.75 40.60 16.72 16.68 5.06 10.93
GBM
duringtheMadonnaproject[48].Thesub-imageofinterest,of “Madonna” images described in paragraph V-A. Then, we
sizeof50×50pixels,isknowntobemainlycomposedofoak fitted a low rank approximation MA to the remaining pixels
and chestnut trees, with an additional unknown non-planted- by minimizing
tree endmember (referred to as Endm. #3 in what follows).
(cid:88)
This dataset will be referred to as the “Madonna” image in d (y |[MA] ) (30)
β lp lp
the following. (l,p)∈O
withrespecttoMandA,whereOdenotesthesetofobserved
B. Selection of β via induction entries. The outlier term R is omitted in this experiment as
Most of the literature in hyperspectral unmixing uses the it cannot be inferred for the missing entries (there is one
squared Euclidean distance for the data-fitting term. This outlier entry per missing data entry and the problem is not
choiceisoftendrivenbycommonpracticeratherthanphysical identifiable). The objective function (30) can be minimized
motivations stemming from the nature of the data. As men- using a minor modification of the MM algorithm described in
tionedinparagraphII-C,divergencesareoftenlog-likelihoods Section III, similarly to the factorizations with missing data
in disguise, and as such, choosing a divergence is akin to described in [34], [49].
making a noise assumption. Thus, when no obvious physical After estimation, the missing pixels (l,p) belonging to
modelsupportsthechoiceofaspecificdivergence,findingthe the complement of O are reconstructed as [MA]lp and the
“right” measure of fit can be seen as a model selection prob- aSAMs values between the original (complete) data Y and its
lem. When data with a ground truth is available for a specific reconstruction Yˆ are computed. This process is repeated for
task,amodelcanbeselectedbasedonitsperformanceforthis values of β from −1 to 3 with a step-size of 0.5. For every
task.Unfortunately,andtothebestofourknowledge,nosuch value of β, 10 runs are considered, corresponding to different
publicrealhyperspectraldataexistsforspectralunmixing,i.e., randominitializationsanddifferentsetsofmissingpixels.The
with perfectly known endmember spectra and corresponding number of endmembers was set to K =3.
abundance coefficients. Another way of selecting a model can The results of the experiment are displayed in Figure 1.
bebasedonitsabilitytopredictunseenormissingdata.Such They show that: (1) the choice of β matters, (2) best per-
an approach does not require a ground truth. As such, in this formance is achieved for β = 1 for the Moffett image and
paragraphwestudyhowNMFwiththeβ-divergenceperforms β =1.5 for the Madonna image, with values of β in the [0,2]
on an interpolation task. Pixels are randomly removed from range yielding sensibly similar performance. The conclusion
real hyperspectral images and those pixels are reconstructed of this study is that it can be worth using alternatives to the
from the low-rank approximation. The process is repeated for standard squared Euclidean distance for hyperspectral unmix-
various values of β and an optimal value of β (in terms of ing (such as the KL divergence, corresponding to β = 1), if
predicting performance) is deduced. it does not come with extra difficulties in the optimization (as
More precisely, 25, 50 or 75% of the pixels have been such, the MM algorithm is equally simple to implement for
randomly and uniformly removed from the ‘Moffett” and all values of β).
8
Vegetation Water Soil
MOFFETT MADONNA
0.4 0.034
p=0.25 p=0.25
0.3
M 0.032
SA0.2
0.1 0.03

−1 0 1 2 3 −1 0 1 2 3
0.4 0.034
p=0.5 p=0.5 0 1 2 0 1 2 0 1 2
0.3
M 0.032
SA0.2
0.1 0.03

−1 0 1 2 3 −1 0 1 2 3
0.4 0.034
p=0.75 p=0.75
0.3
M 0.032
SA0.2 (a) Moffett
0.1 0.03
Chesnut tree Oak tree Endm. #3
−1 0 1 2 3 −1 0 1 2 3
beta beta
Fig. 1. Average SAM values (± standard deviation) between the original
and reconstructed data over 10 runs. The percentage of observed entries is
increased from 25% (top) to 75% (bottom). Left: Moffett image; Right:
Madonnaimage.Bestreconstructionsareobtainedforeitherβ=1or1.5.
0.4 0.6 0.8 1 0.4 0.6 0.8 1 0.4 0.6 0.8 1
C. Robust unmixing results
In a last experiment, the proposed rLMM-unmixing tech-
nique has been applied on the real Moffett and Madonna
images. We used K = 3 and considered two values of β,
namely β = 1 (Kullback-Leibler divergence) and β = 2
(squared Euclidean distance). The endmember spectra and
(b) Madonna
abundancemapsestimatedbyrNMFaredepictedinFig.2.For
conciseness,onlytheabundancemapsobtainedwithβ =1are
Fig.2. Unmixingresultsoftworealhyperspectralimages.Topofeachimage:
displayed as the results for β =2 were visually very similar.
endmembersestimatedbytheproposedrNMF-basedunmixingalgorithmwith
The unmixing results are in good agreement with previ- β=1(redlines)andβ=2(blacklines).Bottomofeachimage:estimated
ous results [47], [50]. However, in addition to the standard abundancemapsobtainedforβ=1;black(resp.white)pixelscorrespondto
absence(resp.presence)oftheassociatedendmembers.
description of the data by linearly mixed endmembers, the
proposedmodelalsoprovidesinformationregardingthepixels
that cannot be explained with the standard LMM. As such,
Fig. 3 displays the energy e = [(cid:107)r (cid:107) ,...,(cid:107)r (cid:107) ] of the
1 2 P 2
residualcomponentestimatedbyrNMF.RegardingtheMoffett
image, the maps demonstrate that most of the pixels of this
scene can be accurately described using the LMM. However,
some few pixels, mainly located in the lake shore, appear
at outliers. These pixels probably correspond to areas where
some interactions between several endmembers occur (e.g.,
water/vegetation,water/soil).Similarresultshavebeenalready
observed in [13], [51], which confirms the relevance of the (a) Moffett (b) Madonna
proposed method. For the Madonna image, the energy map
exhibits residual terms that are mainly located in the area Fig.3. EnergyofthenonlinearcomponentsreturnedbyrNMFwithβ=1.
Dark(resp.light)pixelscorrespondtosmall(resp.large)values.
occupied by the oak trees and the unknown 3rd endmember.
Furthermore,theimageshowsregularverticalpatternsthatare
almost surely due to a sensor defect or miscalibration during
cancaptureso-callednonlineareffects.Thesenonlineareffects
the data post-processing.
are treated as additive and sparsely active outliers. In contrast
with state-of-the-art literature on nonlinear hyperspectral un-
VI. CONCLUSION mixing, our approach does not require the specification of a
In this paper we have presented a new mixing model to particular model of nonlinearity.
describe hyperspectral data. This model, denoted as rLMM, The resulting unmixing problem was formulated as a new
extends the standard LMM by including a residual term that form of robust NMF problem, for which we developed a
9
simple and effective block-coordinate descent algorithm that [16] N. Dobigeon, L. Tits, B. Somers, Y. Altmann, and P. Coppin, “A
involves multiplicative updates. We provided an effective rule comparison of nonlinear mixing models for vegetated areas using
simulated and real hyperspectral data,” IEEE J. Sel. Topics Appl.
of thumb for setting the value of the penalty weight, which
Earth Observations and Remote Sens., submitted. [Online]. Available:
leaves our algorithm virtually free of parameters (only the http://arxiv.org/abs/1312.1270/
number of endmembers needs to be specified). Simulations [17] N.DobigeonandC.Fe´votte,“Robustnonnegativematrixfactorization
fornonlinearunmixingofhyperspectralimages,”inProc.IEEEWork-
conducted on synthetic and real data have illustrated the
shop Hyperspectral image and signal processing: Evolution in remote
relevance of rLMM, which outperformed many unmixing sensing(WHISPERS),Gainesville,FL,2013.
methods designed for various linear and nonlinear models. [18] E.J.Cande`s,X.Li,Y.Ma,andJ.Wright,“Robustprincipalcomponent
analysis?”JournalofACM,vol.58,no.1,pp.1–37,2009.
ACKNOWLEDGEMENTS [19] P.Sprechmann,A.Bronstein,andG.Sapiro,“Real-timeonlinesinging
voiceseparationfrommonauralrecordingsusingrobustlow-rankmod-
We thank Vincent Y. F. Tan and Zhao Renbo from National eling,” in Proc. Int. Soc. Music Information Retrieval Conf. (ISMIR),
UniversityofSingaporefordiscussionsrelatedtorobustNMF Porto,Portugal,Oct.2012.
[20] L.Zhang,Z.Chen,M.Zheng,andX.He,“Robustnonnegativematrix
and feedback about this manuscript.
factorization,” Front. Electr. Electron. Eng. China, vol. 6, no. 2, pp.
192–200,2011.
REFERENCES [21] B. Shen, L. Si, R. Ji, and B. Liu, “Robust nonnegative matrix
[1] G. P. Asner and K. B. Heidebrecht, “Spectral unmixing of vegetation, factorization via (cid:96)1 norm regularization,” ArXiv preprint, 2012.
[Online].Available:http://arxiv.org/abs/1204.2311/
soilanddrycarboncoverinaridregions:comparingmultispectraland
[22] D. Kong, C. Ding, and H. Huang, “Robust nonnegative matrix factor-
hyperspectral observations,” Int. J. Remote Sens., vol. 23, no. 19, pp.
3939–3958,Oct.2002. izationusing(cid:96)21-norm,”inProc.20thACMInt.Conf.Informationand
KnowledgeManagement,2011,pp.673–682.
[2] K. E. Themelis, F. Schmidt, O. Sykioti, A. A. Rontogiannis, K. D.
[23] A.BenHamzaandD.J.Brady,“Reconstructionofreflectancespectra
Koutroumbas, and I. A. Daglis, “On the unmixing of MEx/OMEGA
using robust nonnegative matrix factorizations,” IEEE Trans. Signal
hyperspectral data,” Planetary and Space Science, vol. 68, no. 1, pp.
Process.,vol.54,pp.3637–3642,2006.
34–41,2012.
[3] A.Gowen,C.O’Donnell,P.Cullen,G.Downey,andJ.Frias,“Hyper- [24] V.P.Pauca,J.Piper,andR.J.Plemmons,“Nonnegativematrixfactor-
spectralimaging:anemergingprocessanalyticaltoolforfoodquality izationforspectraldataanalysis,”LinearAlgebraanditsApplications,
and safety control,” Trends in Food Science & Technology, vol. 18, vol.416,pp.29–47,2006.
no.12,pp.590–598,2007. [25] L. Miao and H. Qi, “Endmember extraction from highly mixed data
[4] N.DobigeonandN.Brun,“SpectralmixtureanalysisofEELSspectrum- using minimum volume constrained nonnegative matrix factorization,”
images,”Ultramicroscopy,vol.120,pp.25–34,Sept.2012. IEEETrans.Geosci.andRemoteSensing,vol.45,no.3,pp.765–777,
[5] N. Keshava and J. F. Mustard, “Spectral unmixing,” IEEE Signal 2007.
Process.Mag.,vol.19,no.1,pp.44–57,Jan.2002. [26] Z. Yang, G. Zhou, S. Xie, S. Ding, J.-M. Yang, and J. Zhang, “Blind
[6] J.M.Bioucas-Dias,A.Plaza,N.Dobigeon,M.Parente,Q.Du,P.Gader, spectral unmixing based on sparse nonnegative matrix factorization,”
and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, IEEETrans.ImageProcess.,vol.20,no.4,pp.1112–1125,2011.
statistical,andsparseregression-basedapproaches,”IEEEJ.Sel.Topics [27] E.Esser,M.Moller,S.Osher,G.Sapiro,andJ.Xin,“Aconvexmodel
Appl.EarthObservationsandRemoteSens.,vol.5,no.2,pp.354–379, for nonnegative matrix factorization and dimensionality reduction on
April2012. physicalspace,”IEEETrans.ImageProcess.,vol.21,no.7,pp.3239–
[7] N. Dobigeon, J.-Y. Tourneret, C. Richard, J. C. M. Bermudez, 3252,2012.
S.McLaughlin,andA.O.Hero,“Nonlinearunmixingofhyperspectral [28] B. Jørgensen, “Exponential dispersion models,” J. Roy. Stat. Soc. Ser.
images: Models and algorithms,” IEEE Signal Process. Mag., vol. 31, B,vol.49,no.2,pp.127–162,1987.
no.1,pp.89–94,Jan.2014. [29] A.Basu,I.R.Harris,N.L.Hjort,andM.C.Jones,“Robustandefficient
[8] B.Somers,K.Cools,S.Delalieux,J.Stuckens,D.V.derZande,W.W. estimation by minimising a density power divergence,” Biometrika,
Verstraeten, and P. Coppin, “Nonlinear hyperspectral mixture analysis vol.85,no.3,pp.549–559,Sept.1998.
for tree cover estimates in orchards,” Remote Sens. Environment, vol. [30] A. Cichocki and S. Amari, “Families of Alpha- Beta- and Gamma-
113,pp.1183–1193,Feb.2009. divergences: Flexible and robust measures of similarities,” Entropy,
[9] B. Somers, L. Tits, and P. Coppin, “Quantifying nonlinear spectral vol.12,no.6,pp.1532–1568,June2010.
mixing in vegetated areas: computer simulation model validation and
[31] R.Kompass,“Ageneralizeddivergencemeasurefornonnegativematrix
firstresults,”IEEEJ.Sel.TopicsAppl.EarthObservationsandRemote
factorization,”NeuralComputation,vol.19,no.3,pp.780–791,2007.
Sens.,2014,toappear.
[32] A.Cichocki,R.Zdunek,andS.Amari,“Csiszar’sdivergencesfornon-
[10] Y. Altmann, N. Dobigeon, and J.-Y. Tourneret, “Bilinear models for
negative matrix factorization: Family of new algorithms,” in Proc. 6th
nonlinear unmixing of hyperspectral images,” in Proc. IEEE GRSS
InternationalConferenceonIndependentComponentAnalysisandBlind
Workshop Hyperspectral Image SIgnal Process.: Evolution in Remote
SignalSeparation(ICA),CharlestonSC,USA,Mar.2006,pp.32–39.
Sens.(WHISPERS),Lisbon,Portugal,June2011,pp.1–4.
[33] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono, and
[11] J.M.P.NascimentoandJ.M.Bioucas-Dias,“Nonlinearmixturemodel
S. Sagayama, “Convergence-guaranteed multiplicative algorithms for
forhyperspectralunmixing,”inProc.SPIEImageandSignalProcessing
non-negativematrixfactorizationwithbeta-divergence,”inProc.IEEE
forRemoteSensingXV,L.Bruzzone,C.Notarnicola,andF.Posa,Eds.,
WorkshopMach.LearningforSignalProcess.(MLSP),Sept.2010.
vol.7477,no.1. SPIE,2009,p.74770I.
[34] C.Fe´votteandJ.Idier,“Algorithmsfornonnegativematrixfactorization
[12] W. Fan, B. Hu, J. Miller, and M. Li, “Comparative study between a
withthebeta-divergence,”NeuralComputation,vol.23,no.9,pp.2421–
newnonlinearmodelandcommonlinearmodelforanalysinglaboratory
2456,Sept.2011.
simulated-foresthyperspectraldata,”Int.J.RemoteSens.,vol.30,no.11,
pp.2951–2962,June2009. [35] M. Tweedie, “An index which distinguishes between some important
[13] A. Halimi, Y. Altmann, N. Dobigeon, and J.-Y. Tourneret, “Nonlinear exponential families,” in Proc. Indian Stat. Inst. Golden Jubilee Int.
unmixingofhyperspectralimagesusingageneralizedbilinearmodel,” Conf.,ser.Statistics:ApplicationsandNewDirections,J.K.Ghoshand
IEEE Trans. Geosci. and Remote Sensing, vol. 49, no. 11, pp. 4153– J.Roy,Eds.,Calcutta,India,1984,pp.579–604.
4162,Nov.2011. [36] V. Y. F. Tan and C. Fe´votte, “Automatic relevance determination in
[14] I.Meganem,P.De´liot,X.Briottet,Y.Deville,andS.Hosseini,“Linear- nonnegativematrixfactorizationwiththebeta-divergence,”IEEETrans.
quadratic mixing model for reflectances in urban environments,” IEEE Patt.Anal.Mach.Intell.,vol.35,no.7,pp.1592–1605,July2013.
Trans. Geosci. and Remote Sensing, vol. 52, no. 1, pp. 544–558, Jan. [37] C. Fe´votte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix factor-
2014. ization with the Itakura-Saito divergence. With application to music
[15] Y.Altmann,A.Halimi,N.Dobigeon,andJ.-Y.Tourneret,“Supervised analysis,”NeuralComputation,vol.21,no.3,pp.793–830,Mar.2009.
nonlinear spectral unmixing using a post-nonlinear mixing model for [38] J. Eggert and E. Ko¨rner, “Sparse coding and NMF,” in Proc. IEEE
hyperspectralimagery,”IEEETrans.ImageProcess.,vol.21,no.6,pp. International Joint Conference on Neural Networks, 2004, pp. 2529–
3017–3025,June2012. 2533.
10
[39] T.Virtanen,“Monauralsoundsourceseparationbynon-negativematrix NicolasDobigeonreceivedthestateengineeringdegreeinelectricalengineer-
factorization with temporal continuity and sparseness criteria,” IEEE ingfromENSEEIHT,Toulouse,France,andtheM.Sc.degreeinsignalpro-
Trans.Audio,Speech,LanguageProcess.,vol.15,no.3,pp.1066–1074, cessingfromtheNationalPolytechnicInstituteofToulouse(INPToulouse),
Mar.2007. bothinJune2004,aswellasthePh.D.degreeandHabilitationa`Dirigerdes
[40] A.Lee,F.Caron,A.Doucet,andC.Holmes,“AHierarchicalBayesian Recherches in Signal Processing from the INP Toulouse in 2007 and 2012,
FrameworkforConstructingSparsity-inducingPriors,”arXiv.org,Sept. respectively.HewasaPost-DoctoralResearchAssociatewiththeDepartment
2010. of Electrical Engineering and Computer Science, University of Michigan,
[41] RSI(ResearchSystemsInc.),ENVIUser’sguideVersion4.0,Boulder, Ann Arbor, MI, USA, from 2007 to 2008. Since 2008, he has been with
CO80301USA,Sept.2003. the National Polytechnic Institute of Toulouse (INP-ENSEEIHT, University
[42] J.M.NascimentoandJ.M.Bioucas-Dias,“Vertexcomponentanalysis: of Toulouse) where he is currently an Associate Professor. He conducts his
afastalgorithmtounmixhyperspectraldata,”IEEETrans.Geosci.and researchwithintheSignalandCommunicationsGroupoftheIRITLaboratory
RemoteSensing,vol.43,no.4,pp.898–910,April2005. and he is also an affiliated faculty member of the Telecommunications for
[43] D. C. Heinz and C. -I Chang, “Fully constrained least-squares linear Space and Aeronautics (TeSA) cooperative laboratory. His current research
spectral mixture analysis method for material quantification in hyper- interests include statistical signal and image processing, with a particular
spectral imagery,” IEEE Trans. Geosci. and Remote Sensing, vol. 29, interest in Bayesian inverse problems with applications to remote sensing,
no.3,pp.529–545,March2001. biomedicalimagingandgenomics.
[44] R. Heylen, D. Burazerovic, and P. Scheunders, “Non-linear spectral
unmixing by geodesic simplex volume maximization,” IEEE J. Sel.
TopicsSignalProcess.,vol.5,no.3,pp.534–542,June2011.
[45] A. Halimi, Y. Altmann, N. Dobigeon, and J.-Y. Tourneret, “Unmixing
hyperspectral images using the generalized bilinear model,” in Proc.
IEEE Int. Conf. Geosci. Remote Sens. (IGARSS), Vancouver, Canada,
July2011,pp.1886–1889.
[46] JetPropulsionLab.(JPL),“Avirisfreedata,”CaliforniaInst.Technol.,
Pasadena,CA,2006.[Online].Available:http://aviris.jpl.nasa.gov/html/
aviris.freedata.html
[47] N.Dobigeon,J.-Y.Tourneret,andC.-IChang,“Semi-supervisedlinear
spectralunmixingusingahierarchicalBayesianmodelforhyperspectral
imagery,”IEEETrans.SignalProcess.,vol.56,no.7,pp.2684–2695,
July2008.
[48] D.Sheeren,M.Fauvel,S.Ladet,A.Jacquin,G.Bertoni,andA.Gibon,
“Mapping ash tree colonization in an agricultural mountain landscape:
Investigatingthepotentialofhyperspectralimagery,”inProc.IEEEInt.
Conf. Geosci. Remote Sens. (IGARSS), Vancouver, Canada, July 2011,
pp.3672–3675.
[49] N.-D. Ho, “Nonnegative matrix factorization algorithms and
applications,” Ph.D. dissertation, Universit Catholique de Louvain,
2008.[Online].Available:www.inma.ucl.ac.be/∼vdooren/ThesisHo.pdf
[50] Y. Altmann, N. Dobigeon, S. McLaughlin, and J.-Y. Tourneret, “Non-
linear spectral unmixing of hyperspectral images using Gaussian pro-
cesses,” IEEE Trans. Signal Process., vol. 61, no. 10, pp. 2442–2453,
May2013.
[51] O.Besson,N.Dobigeon,andJ.-Y.Tourneret,“Minimummeansquare
distanceestimationofasubspace,”IEEETrans.SignalProcess.,vol.59,
no.12,pp.5709–5720,Dec.2011.
Ce´dric Fe´votte received the state engineering and PhD degrees in control
and computer science from the E´cole Centrale de Nantes, France, in 2000
and 2003, respectively. During his PhD, he was with the Signal Processing
Group at the Institut de Recherche en Communication et Cyberne´tique de
Nantes(IRCCyN).From2003to2006,hewasaresearchassociatewiththe
Signal Processing Laboratory at the University of Cambridge (Engineering
Department).Hewasthenaresearchengineerwiththemusiceditingtechnol-
ogystart-upcompanyMist-Technologies(nowAudionamix)inParis.In2007,
hebecameaCNRStenuredresearcher.HewasaffiliatedwithLTCI(CNRS
& Te´le´com ParisTech) from 2007 to 2012. Since 2013, he has been with
LaboratoireLagrangre(CNRS,ObservatoiredelaCoˆted’Azur&Universite´
deNiceSophiaAntipolis).Hisresearchinterestsgenerallyconcernstatistical
signalprocessingandmachinelearning,inparticularforinverseproblemsand
sourceseparation.HeisamemberoftheIEEE“MachineLearningforSignal
Processing”technicalcommittee.