2014 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 49, NO. 6, JUNE 2011
Sparse Unmixing of Hyperspectral Data
Marian-Daniel Iordache, José M. Bioucas-Dias, Member, IEEE, and Antonio Plaza, Senior Member, IEEE
Abstract—Linear spectral unmixing is a popular tool in re-
motely sensed hyperspectral data interpretation. It aims at esti-
mating the fractional abundances of pure spectral signatures (also
called as endmembers) in each mixed pixel collected by an imaging
spectrometer. In many situations, the identification of the end-
member signatures in the original data set may be challenging due
to insufficient spatial resolution, mixtures happening at different
scales, and unavailability of completely pure spectral signatures in
the scene. However, the unmixing problem can also be approached
in semisupervised fashion, i.e., by assuming that the observed im-
age signatures can be expressed in the form of linear combinations
of a number of pure spectral signatures known in advance (e.g.,
spectra collected on the ground by a field spectroradiometer).
Unmixing then amounts to finding the optimal subset of signatures
in a (potentially very large) spectral library that can best model
each mixed pixel in the scene. In practice, this is a combina-
torial problem which calls for efficient linear sparse regression
(SR) techniques based on sparsity-inducing regularizers, since the
number of endmembers participating in a mixed pixel is usually
very small compared with the (ever-growing) dimensionality (and
availability) of spectral libraries. Linear SR is an area of very
active research, with strong links to compressed sensing, basis
pursuit (BP), BP denoising, and matching pursuit. In this paper,
we study the linear spectral unmixing problem under the light
of recent theoretical results published in those referred to areas.
Furthermore, we provide a comparison of several available and
new linear SR algorithms, with the ultimate goal of analyzing their
potential in solving the spectral unmixing problem by resorting to
available spectral libraries. Our experimental results, conducted
using both simulated and real hyperspectral data sets collected by
the NASA Jet Propulsion Laboratory’s Airborne Visible Infrared
Imaging Spectrometer and spectral libraries publicly available
from the U.S. Geological Survey, indicate the potential of SR
techniques in the task of accurately characterizing the mixed
pixels using the library spectra. This opens new perspectives for
spectral unmixing, since the abundance estimation process no
longer depends on the availability of pure spectral signatures in the
input data nor on the capacity of a certain endmember extraction
algorithm to identify such pure signatures.
Index Terms—Abundance estimation, convex optimization,
hyperspectral imaging, sparse regression (SR), spectral unmixing.
Manuscript received March 14, 2010; revised August 18, 2010 and
November 2, 2010; accepted November 28, 2010. Date of publication
January 19, 2011; date of current version May 20, 2011. This work was
supported by the European Community’s Marie Curie Research Training
Networks Program under contract MRTN-CT-2006-035927 (Hyperspectral
Imaging Network) and by the Spanish Ministry of Science and Innovation
(HYPERCOMP/EODIX project, reference AYA2008-05965-C04-02).
M.-D. Iordache and A. Plaza are with the Department of Technology of Com-
puters and Communications, Escuela Politécnica, University of Extremadura,
10071 Cáceres, Spain (e-mail: diordache@unex.es; aplaza@unex.es).
J. M. Bioucas-Dias is with the Instituto de Telecomunicações, Instituto,
Superior Técnico, 1049-1 Lisbon, Portugal (e-mail: bioucas@lx.it.pt).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TGRS.2010.2098413
I. INTRODUCTION
H
YPERSPECTRAL imaging has been transformed from
being a sparse research tool into a commodity product
that is available to a broad user community [1]. The wealth
of spectral information available from advanced hyperspectral
imaging instruments currently in operation has opened new
perspectives in many application domains, such as monitoring
of environmental and urban processes or risk prevention and
response, including, among others, tracking wildfires, detecting
biological threats, and monitoring oil spills and other types of
chemical contamination. Advanced hyperspectral instruments
such as NASA’s Airborne Visible Infrared Imaging Spectrom-
eter (AVIRIS) [2] are now able to cover the wavelength region
from 0.4 to 2.5 μm using more than 200 spectral channels at
a nominal spectral resolution of 10 nm. The resulting hyper-
spectral data cube is a stack of images (see Fig. 1) in which
each pixel (vector) is represented by a spectral signature or
fingerprint that characterizes the underlying objects.
Several analytical tools have been developed for remotely
sensed hyperspectral data processing in recent years, cover-
ing topics like dimensionality reduction, classification, data
compression, or spectral unmixing [3], [4]. The underlying
assumption governing clustering and classification techniques
is that each pixel vector comprises the response of a single
underlying material. However, if the spatial resolution of the
sensor is not high enough to separate different materials, these
can jointly occupy a single pixel. For instance, it is likely that
the pixel collected over a vegetation area in Fig. 1 actually
comprises a mixture of vegetation and soil. In this case, the
measured spectrum may be decomposed into a linear combina-
tion of pure spectral signatures of soil and vegetation, weighted
by abundance fractions that indicate the proportion of each
macroscopically pure signature in the mixed pixel [5].
To deal with this problem, linear spectral mixture analysis
techniques first identify a collection of spectrally pure con-
stituent spectra, called as endmembers in the literature, and then
express the measured spectrum of each mixed pixel as a linear
combination of endmembers weighted by fractions or abun-
dances that indicate the proportion of each endmember present
in the pixel [6]. It should be noted that the linear mixture model
assumes minimal secondary reflections and/or multiple scatter-
ing effects in the data collection procedure, and hence, the mea-
sured spectra can be expressed as a linear combination of the
spectral signatures of the materials present in the mixed pixel
[see Fig. 2(a)]. Being quite opposite, the nonlinear mixture
model assumes that the endmembers form an intimate mixture
inside the respective pixel so that the incident radiation interacts
with more than one component and is affected by multiple
scattering effects [see Fig. 2(b)]. Nonlinear unmixing generally
0196-2892/$26.00 © 2011 IEEE
IORDACHE et al.: SPARSE UNMIXING OF HYPERSPECTRAL DATA 2015
Fig. 1. Concept of hyperspectral imaging and the presence of mixed pixels.
Fig. 2. (a) Linear versus (b) nonlinear mixture models.
requires prior knowledge about object geometry and the phys-
ical properties of the observed objects. In this paper, we will
focus exclusively on the linear mixture model due to its com-
putational tractability and flexibility in different applications.
The linear mixture model assumes that the spectral response
of a pixel in any given spectral band is a linear combination
of all of the endmembers present in the pixel at the respective
spectral band. For each pixel, the linear model can be written as
follows:
y
i
=
q
j=1
m
ij
α
j
+ n
i
(1)
where y
i
is the measured value of the reflectance at spectral
band i, m
ij
is the reflectance of the jth endmember at spectral
band i, α
j
is the fractional abundance of the jth endmember,
and n
i
represents the error term for the spectral band i (i.e., the
noise affecting the measurement process). If we assume that
the hyperspectral sensor used in data acquisition has L spectral
bands, (1) can be rewritten in compact matrix form as
y = Mα + n (2)
where y is an L × 1 column vector (the measured spectrum
of the pixel), M is an L × q matrix containing q pure spec-
tral signatures (endmembers), α is a q × 1 vector containing
the fractional abundances of the endmembers, and n is an
L × 1 vector collecting the errors affecting the measurements
at each spectral band. The so-called abundance nonnegativity
constraint (ANC) (α
i
≥ 0 for i =1,...,q) and the abundance
sum-to-one constraint (ASC) (
q
i=1
α
i
=1), which we, re-
spectively, represent in compact form by
α ≥0 (3)
1
T
α =1 (4)
where 1
T
is a line vector of 1’s compatible with α,areoften
imposed into the model described in (1) [7], owing to the
fact that α
i
,fori =1,...,q, represents the fractions of the
endmembers present in the considered pixel.
In a typical hyperspectral unmixing scenario, we are given a
set Y ≡{y
i
∈ R
L
,i=1,...,n}of n observed L-dimensional
spectral vectors, and the objective is to estimate the mixing
matrix M and the fractional abundances α for every pixel in the
scene. This is a blind source separation problem, and naturally,
independent component analysis methods come to mind to
solve it. However, the assumption of statistical independence
among the sources (the fractional abundances in our applica-
tion), central to independent component analysis methods, does
not hold in hyperspectral applications, since the sum of frac-
tional abundances associated to each pixel is constant. Thus,
the sources are statistically dependent, which compromises the
performance of independent component analysis algorithms in
hyperspectral unmixing [8].
2016 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 49, NO. 6, JUNE 2011
We note that the constraints (3) and (4) define the set S
q−1
≡
{α ∈ R
q
|α ≥ 0, 1
T
α =1}, which is the probability simplex
in R
q
. Furthermore, the set S
M
≡{Mα ∈ R
L
|α ∈S
q−1
} is
also a simplex whose vertices are the columns of M.Overthe
last decade, several algorithms have exploited this geometrical
property by estimating the “smallest” simplex set containing
the observed spectral vectors [9], [10]. Some classic techniques
for this purpose assume that the input data set contains at
least one pure pixel for each distinct material present in the
scene, and therefore, a search procedure aimed at finding the
most spectrally pure signatures in the input scene is feasible.
Among the endmember extraction algorithms working under
this regime, we can list some popular approaches such as the
pixel purity index [11], N-FINDR [12], orthogonal subspace
projection technique in [13], and vertex component analysis
(VCA) [14]. However, the assumption under which these al-
gorithms perform may be difficult to guarantee in practical
applications due to several reasons.
1) First, if the spatial resolution of the sensor is not high
enough to separate different pure signature classes at a
macroscopic level, the resulting spectral measurement
can be a composite of individual pure spectra which cor-
respond to materials that jointly occupy a single pixel. In
this case, the use of image-derived endmembers may not
result in accurate fractional abundance estimations since
it is likely that such endmembers may not be completely
pure in nature.
2) Second, mixed pixels can also result when distinct mate-
rials are combined into a microscopic (intimate) mixture,
independent from the spatial resolution of the sensor.
Since the mixtures in this situation happen at the parti-
cle level, the use of image-derived spectral endmembers
cannot accurately characterize intimate spectral mixtures.
In order to overcome the two aforementioned issues, other
advanced endmember generation algorithms have also been
proposed under the assumption that pure signatures are not
present in the input data. Such techniques include optical
real-time adaptive spectral identification systems [15], convex
cone analysis [16], iterative error analysis [17], automatic
morphological endmember extraction [18], iterated constrained
endmembers (ICE) [19], minimum volume constrained non-
negative matrix factorization [20], spatial–spectral endmember
extraction [21], sparsity-promoting ICE [22], minimum vol-
ume simplex analysis [23], and simplex identification via split
augmented Lagrangian [24]. A necessary condition for these
endmember generation techniques to yield good estimates is the
presence in the data set of at least q − 1 spectral vectors on each
facet of the simplex set S
M
[24]. This condition is very likely
to fail in highly mixed scenarios, in which the aforementioned
techniques generate artificial endmembers, i.e., not necessarily
associated to physically meaningful spectral signatures of true
materials.
In this paper, we adopt a novel semisupervised approach
to linear spectral unmixing, which relies on the increasing
availability of spectral libraries of materials measured on the
ground, e.g., using advanced field spectroradiometers. Our
main assumption is that mixed pixels can be expressed in
the form of linear combinations of a number of pure spectral
signatures known in advance and available in a library, such as a
the well-known one publicly available from the U.S. Geological
Survey (USGS),
1
which contains over 1300 mineral signatures,
or the NASA Jet Propulstion Laboratory’s Advanced Space-
borne Thermal Emission and Reflection Radiometer (ASTER)
spectral library,
2
which is a compilation of over 2400 spectra of
natural and man-made materials. When the unmixing problem
is approached using spectral libraries, the abundance estimation
process no longer depends on the availability of pure spectral
signatures in the input data nor on the capacity of a certain
endmember extraction algorithm to identify such pure signa-
tures. Being quite opposite, the procedure is reduced to finding
the optimal subset of signatures in the library that can best
model each mixed pixel in the scene. Despite the appeal of this
semisupervised approach to spectral unmixing, this approach is
also subject to a few potential drawbacks.
1) One risk in using library endmembers is that these spectra
are rarely acquired under the same conditions as the
airborne data. Image endmembers have the advantage of
being collected at the same scale as the data, and thus,
they can be more easily associated with features on the
scene. However, such image endmembers may not always
be present in the input data. In this paper, we rely on
the use of advanced atmospheric correction algorithms
which convert the input hyperspectral data from at-sensor
radiance to reflectance units.
2) The ability to obtain useful sparse solutions for an un-
derdetermined system of equations mostly depends on
the degree of coherence between the columns of the
system matrix and the degree of sparseness of the original
signals (i.e., the abundance fractions) [25]–[28]. The most
favorable scenarios correspond to highly sparse signals
and system matrices with low coherence. Unfortunately,
in hyperspectral applications, the spectral signatures of
the materials tend to be highly correlated. On the other
hand, the number of materials present in a given scene
is often small, e.g., less than 20, and most importantly,
the number of materials participating in a mixed pixel
is usually on the order of four to five [5]. Therefore, the
undesirable high coherence of hyperspectral libraries can
be mitigated, to some extent, by the highly sparse nature
of the original signals.
3) The sparse solutions of the underdetermined systems are
computed by solving the optimization problems contain-
ing nonsmooth terms [26]. The presence of these terms
introduces complexity because the standard optimization
tools of the gradient and Newton family cannot be di-
rectly used. To make the scenario even more complex,
a typical hyperspectral image has hundreds or thousands
of spectral vectors, implying an equal number of inde-
pendent optimizations to unmix the complete scene. To
cope up with this computational complexity, we resort
to recently introduced (fast) algorithms based on the
augmented Lagrangian method of multipliers [29].
1
Available online at http://speclab.cr.usgs.gov/spectral-lib.html.
2
Available online at http://speclib.jpl.nasa.gov.
IORDACHE et al.: SPARSE UNMIXING OF HYPERSPECTRAL DATA 2017
In this paper, we specifically address the problem of spar-
sity when unmixing the hyperspectral data sets using spectral
libraries and further provide a quantitative and comparative as-
sessment of several available and new optimization algorithms
in the context of linear sparse problems. The remainder of this
paper is organized as follows. Section II formulates the sparse
regression (SR) problem in the context of hyperspectral unmix-
ing. Section III describes several available and new unmixing
algorithms, with the ultimate goal of analyzing their potential in
solving the sparse hyperspectral unmixing problems. Section IV
provides an experimental validation of the considered algo-
rithms using the simulated hyperspectral mixtures from the real
and synthetic spectral libraries. The primary reason for the use
of the simulated data is that all details of the simulated mixtures
are known, and they can be efficiently investigated because they
can be manipulated individually and precisely. As a comple-
ment to the simulated data experiments, Section V presents an
experimental validation of the considered SR and convex op-
timization algorithms using a well-known hyperspectral scene
collected by the AVIRIS instrument over the Cuprite mining
district in NV. The USGS spectral library is used in conducting
extensive semisupervised unmixing experiments on this scene.
Finally, Section VI concludes with some remarks and hints
at plausible future research. The Appendix is devoted to the
description of the parameter settings used in our experiments
and to the strategies followed to infer these parameters.
II. S
PECTRAL UNMIXING REFORMULATED
AS AN
SR PROBLEM
In this section, we revisit the classic linear spectral unmixing
problem and reformulate it as a semisupervised approach using
SR terminology. Furthermore, we review the SR optimization
problems that are relevant to our unmixing problem, their theo-
retical characterization, their computational complexity, and the
algorithms that are used to solve them exactly or approximately.
Let us assume that the spectral endmembers that are used to
solve the mixture problem are no longer extracted nor generated
using the original hyperspectral data as input but are, instead,
selected from a library containing a large number of spectral
samples available apriori. In this case, unmixing amounts to
finding the optimal subset of samples in the library that can
best model each mixed pixel in the scene. This means that a
searching operation must be conducted in a (potentially very
large) library, which we denote by A ∈ R
L×m
, where L and
m are the number of spectral bands and the number of mate-
rials in the library, respectively. All libraries herein considered
correspond to underdetermined systems, i.e., L<m. With the
aforementioned assumptions in mind, let x ∈ R
m
denote the
fractional abundance vector with regard to the library A.As
usual, we say that x is a k-sparse vector if it has, at most, k
components different from zero. With these definitions in place,
we can now write our SR problem as
min
x
x
0
subject to y−Ax
2
≤δ, x≥0, 1
T
x=1
(5)
where x
0
denotes the number of nonzero components of x
and δ ≥ 0 is the error tolerance due to the noise and modeling
errors. The solution of problem (5), if any, belongs to the set of
sparsest signals belonging to the (m − 1)-probability simplex
satisfying error tolerance inequality y − Ax
2
≤ δ. Prior to
addressing problem (5), we consider a series of simpler related
problems.
A. Exact Solutions
Let us first start by assuming that the noise is zero and the
ANC and ASC constraints are not enforced. Our SR optimiza-
tion problem is then
(P
0
): min
x
x
0
subject to Ax = y. (6)
If the system of linear equations Ax = y has a solution
satisfying 2x
0
< spark(A), where spark(A) ≤ rank(A)+1
is the smallest number of linearly dependent columns of A,itis
necessarily the unique solution of (P
0
) [30], [31]. The spark of
a matrix gives us a very simple way to check the uniqueness of a
solution of the system Ax = y. For example, if the elements of
A are independent and identically distributed (i.i.d.), then with
a probability of one, we have spark(A)=m +1, implying that
every solution with no more than L/2 entries is unique.
In our SR problem, we would like then to compute the spark
of the hyperspectral library being used to have an idea of what
is the minimum level of sparsity of the fractional abundance
vectors that can be uniquely determined by solving (P
0
).Com-
puting the spark of a general matrix is, however, a hard prob-
lem, at least as difficult as solving (P
0
). This complexity has
fostered the introduction of entities that are simpler to compute,
although providing less tight bounds. Mutual coherence is such
an example. Denoting the kth column in A by a
k
and the
2
norm by ·
2
, the mutual coherence of A is given by
μ(A) ≡ max
1≤k,j≤m,k=j
a
T
k
a
j
a
k
2
a
j
2
(7)
i.e., by the maximum absolute value of the cosine of the angle
between any two columns of A. Mutual coherence supplies us
with a lower bound for the spark given by [30]
spark(A) ≥ 1+
1
μ(A)
.
Unfortunately, as it will be shown further, the mutual coherence
of the hyperspectral libraries is very close to one, leading to
useless bounds for the spark. In the following, we illustrate two
relaxed strategies for computing (P
0
): pursuit algorithms and
nonnegative signals.
1) Pursuit Algorithms: The problem (P
0
) is NP hard (which
means that the problem is combinatorial and very complex
to solve) [32], and therefore, there is a little hope in solving
it in a straightforward way. Greedy algorithms such as the
orthogonal basis pursuit [orthogonal matching pursuit (OMP)]
[33] and basis pursuit (BP) [34] are two alternative approaches
in computing the sparsest solution. BP replaces the
0
norm in
(P
0
) with the
1
norm
(P
1
): min
x
x
1
subject to Ax = y. (8)
2018 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 49, NO. 6, JUNE 2011
Contrary to problem (P
0
), problem (P
1
) is convex, and it can
be written as a linear programming (LP) problem and can be
solved using LP solvers. What is, perhaps, totally unexpected is
that, in given circumstances related to matrix A, problem (P
1
)
has the same solution as problem (P
0
). This result is stated in
terms of the restricted isometric constants introduced in [27].
Herein, we use the variant proposed in [35]. Let α
k
, with β
k
≥
0, be the tightest constants in the inequalities
α
k
x
2
≤Ax
2
≤ β
k
x
2
, x
0
≤ k (9)
and further define
γ
2s
≡
β
2
2s
α
2
2s
≥ 1. (10)
Then, under the assumption that γ
2s
< 4
√
2 − 3 2.6569,
every s-sparse vector is recovered by solving problem (P
1
)
[35, Th. 2.1 and Corol. 2.1]. Meanwhile, it has been shown
that, in some cases, the OMP algorithm also provides the (P
0
)
solution in a fashion that is comparable with the BP alternative,
with the advantage of being faster and easier to implement
[26], [36].
2) Nonnegative Signals: We now consider the problem
P
+
0
:min
x
x
0
subject to Ax = yx≥ 0 (11)
and follow a line of reasoning that is close to that of [25]. The
hyperspectral libraries generally contain only the nonnegative
components (i.e., reflectances). Thus, by assuming that the zero
vector is not in the columns of A, it is always possible to find a
vector h such that
h
T
A = w
T
> 0. (12)
Since all of the components of w are nonnegative, matrix W
−1
,
where W ≡ diag(w), is well defined, and it has positive diago-
nal entries. Defining z ≡ Wx, c ≡ h
T
y, and D ≡ AW
−1
and
noting that
h
T
AW
−1
z = 1
T
z (13)
the problem (P
+
0
) is equivalent to
P
+
0
:min
x
z
0
subject to Dz=yz≥0, 1
T
z=c.
(14)
We conclude that, when the original signals are nonnegative and
the system matrices comply with property (12), then problem
(11) enforces the equality constraint 1
T
z = c. This constraint
has very strong connections with the ASC constraint which is
so popular in hyperspectral applications. The ASC is, however,
prone to strong criticisms because, in a real image, there is a
strong signature variability [37] that, at the very least, intro-
duces positive scaling factors varying from pixel to pixel in the
signatures present in the mixtures. As a result, the signatures
are defined up to a scale factor, and thus, the ASC should be
replaced with a generalized ASC of the form
i
ξ
i
x
i
=1,in
which the weights ξ
i
denote the pixel-dependent scale factors.
What we conclude from the equivalence between problems (11)
and (14) is that the nonnegativity of the sources automatically
imposes a generalized ASC. For this reason, we do not explic-
itly impose the ASC constraint.
Similar to problem (P
0
), problem (P
+
0
) is NP hard and
impossible to exactly solve for a general matrix A.Asin
Section II-A1, we can consider instead
1
relaxation
P
+
1
:min
x
z
1
subject to Dz = yz≥ 0. (15)
Here, we have dropped the equality constraint 1
T
z = c because
it is satisfied by any solution of Dz = y. As with problem
(P
0
), the condition γ
2s
< 4
√
2 − 3 2.6569 referred to in
Section II-A1 is now applied to the restricted isometric con-
stants of matrix D to ensure that any s-sparse vector solution of
(P
+
0
) is recovered by solving the problem (P
+
1
).
Another way of characterizing the uniqueness of the solution
of problem (P
+
0
) is via the one-sided coherence introduced
in [25]. However, similar to mutual coherence, the one-sided
coherence of the hyperspectral libraries is very close to one,
leading to useless bounds. The coherence may be increased by
left multiplying the system Dz = y with a suitable invertible
matrix P[25]. This preconditioning tends to improve the per-
formance of greedy algorithms such as OMP. It leads, however,
to an optimization problem that is equivalent to (P
+
1
).Thus,a
BP solver yields the same solution.
B. Approximate Solutions
We now assume that the perturbation n in the observation
model is not zero, and we still want to find an approximate
solution for our SR problem. The computation of the approx-
imate solutions raises issues that are parallel to those found
for exact solutions as addressed earlier. Therefore, we go very
briefly through the same topics. Again, we start by assuming
that the noise is zero and the ANC and ASC constraints are not
enforced. Our noise-tolerant SR optimization problem is then
P
δ
0
:min
x
x
0
subject to Ax − y
2
≤ δ. (16)
The concept of uniqueness of the sparsest solution is now
replaced with the concept of stability [35], [38], [39]. For
example, in [38], it is shown that, given a sparse vector x
0
satisfying the sparsity constraint x
0
< (1 + 1/μ(A))/2 such
that Ax
0
− y≤δ, then every solution x
δ
0
of problem (P
δ
0
)
satisfies
x
δ
0
− x
0
2
≤
4δ
2
1 − μ(A)(2x
0
− 1)
. (17)
Notice that, when δ =0, i.e., when the solutions are exact, this
result parallels those ensuring the uniqueness of the sparsest
solution. Again, we illustrate two relaxed strategies for com-
puting (P
0
).
1) Pursuit Algorithms: Problem (P
δ
0
),as(P
0
), is NP hard.
We consider here two approaches to tackle this problem. The
first approach is the greedy OMP algorithm with stopping rule
Ax − y
2
≤ δ. The second one consists of relaxing the
0
