PARALLEL SPARSE UNMIXING OF HYPERSPECTRAL DATA Jose M. Rodriguez A lvesa, Jose M.
P. Nascimentoa, b , Jose M. Bioucas-Diasa,c, Antonio Plaza d, Vftor Silvae
alnstituto de Teleco munica90es, Lisbon, Portugal b lnstituto Superior de Engenharia de Lisboa, Lisbon, Portugal clnstituto Superior Tecnico , Technical University of Lisbon, Lisbon, Portugal dHyperspectral Computing Laboratory, University of Extremadura, Caceres, Spain elnstituto de Telecomunica90es, DEEC, University of Coimbra, Coimbra, Portugal
ABSTRACT
spectrally distinct materials (also called
In this paper, a new parallel method for sparse spectral un mixing of remotely sensed hyperspectral data on commodity graphics processing units (GPUs) is presented. A semi-supervised approach is adopted, which relies on the increasing availability of spectral libraries of materials measured on the ground instead of resorting to endmember extraction methods. This method is based on the
mixing by splitting and augmented Lagrangian
spectral un
(SUNSAL)
that estimates the material' s abundance fractions.
The par
allel method is performed in a pixel-by-pixel fashion and its implementation properly exploits the GPU architecture at low level, thus taking full advantage of the computational power
endmembers)
can be
found within the same pixel. Thus, each pixel can be viewed as a mixture of endmember signatures . Hyperspectral unmixing is a source separation problem which aims at estimating the number of endmembers, their spectral signatures and their abundance fractions percentage of each endmember)
[1].
(i. e. ,
the
Over the last decade,
several algorithms have been developed to unmix remotely sensed hyperspectral data, from a geometrical and statisti cal point of view
[1].
Some of this methods
[5, 6, 7, 8]
as
sume that dataset contains at least one pure pixel of each dis tinct endmember. However, when the remotely sensed data is highly mixed this assumption may not be valid
[9] .
Experimental results obtained for simulated and
A semi-supervised strategy to cope the unavailability of
real hyperspectral datasets reveal significant speedup factors,
pure spectral signatures is to model mixed pixel observations
of GPUs. up to
1 64
times, with regards to optimized serial implementa
as linear combinations of spectra from a library collected on the ground by a field spectro-radiometer
tion. Index Terms- Hyperspectral Unmixing, Sparse Regres
sion , Graphics Processing Unit, Parallel Methods.
[ 1 0, 1 1 ] .
Thus, un
mixing problem consists in findind the optimal subset of sig natures in a potentially very large spectral library that can best model each mixed pixel in the scene.
In practice, this is a
combinatorial problem which calls for efficient linear sparse
1. INTRODUCTION
regression techniques based on sparsity-inducing regulariz
Remotely sensed hyperspectral images collect electromag
pixel is usually very small compared with the dimensionality
netic energy scattered within their ground instantaneous field
and availability of spectral libraries
of view in hundreds of nearly contiguous spectral bands with
sparse techniques applied to hyperspectral unmixing can be
ers, since the number of endmembers participating in a mixed
high spectral resolution
[1].
This technology provides enough
spectral resolution for material identification, facilitating an enormous number of applications in the fields of urban and regional planning,
water resource management,
environ
mental monitoring, food safety, counterfeit drugs detection, wild land fire tracking, oil spill and other types of chemical contamination detection, biological hazards prevention, and target detection for military and security purposes, among many others
[2, 3 ] .
Due t o low spatial resolution provided b y these devices and to other effects (see This
work
was
[ 1 , 4]
supported
by
for more details), several FeT-IT
under
project
PEst
OE/EEIILA0008120 1 3 .
978-1-4799-1114-1/13/$31.00 ©201 3 IEEE
1446
found in
[ 1 0] .
Some examples of
[ 1 2, 1 3 , 1 4, 1 5 , 1 6] .
However, due to the high dimensionality o f the hyper spectral scene and to the potentially very large library, the aforementioned techniques may be computationally very ex pensive, which prevents their use in time-critical applications
[ 1 7, 1 8] .
Despite the growing interest in parallel hyper
spectral imaging methods and in their GPU implementations
[ 1 9] ,
no parallel implementations of sparse unmixing tech
niques are available on the open literature so far. This paper proposes a parallel method designed for graphics process ing units (GPUs), to solve the constrained sparse regression problem. This method is based on the
splitting and augmented Lagrangian
spectral unmixing by
(SUNSAL)
[20]
that
IGARSS 2013
estimates the abundance fractions using the alternating di rection method of multipliers (ADMM)
[2 1 ] .
The technique
S . This problem is a particular case of the con £2 - £ 1 problems solved by SUNSAL, corresponding absence of the £ 1 term. SUNSAL algorithm uses the
sparsity of strained
decomposes a difficult problem into a sequence of simpler
to the
ones, alleviating considerably the computational burden, and
alternating direction method of multipliers (ADMM)
allowing a significant increase in processing speed.
solve
Addi
(2).
[2 1 ]
to
The pseudo-code is presented in Algorithm 1 .
tionally, this method is performed in a parallel pixel-by-pixel fashion and implemented on a GPU exploiting its architecture at low level, thus taking full advantage of its computational power, extremely high floating-point processing performance, and huge memory bandwidth. The remainder of the paper is organized as follows . Sec tion
2
formulates the problem and presents the fundamentals
of the sparse regression problem. Section 3 describes the pro posed parallel method designed for GPU. Section
4
evaluates
the acceleration of the proposed method from the computa tional point of view. Finally, section 5 concludes the paper with some remarks.
Algorithm 1 2: 3: 4: 5: 6:
repeat
7:
R := H + p, (V k + D k ) S k+ l := GR + cl� V k := Sk+l - D k V k+l := ma.x{O , Vd D k+l := D k - ( S k+ l - V k+l ) k := k + 1
8: 9: 10: 11:
2. SPARSE UNMIXING FORMULATION
12:
Linear mixing model considers that a I - dimensional mixed pixel
y,
of an hyperspectral image with I spectral bands, is a
SUNSAL
:
p, > 0, Va, and Do. H := ATy B := A T A + p,I 1 c := B - l 1m (1;z:,B- l 1mr l 1 G := B - - c1;;; B k := °
choose
I:
13:
until stopping criterion is satisfied
14:
linear combination of endmember signatures weighted by the correspondent abundance fractions.
m
Let A
==
[a l , a2 , . . . , am]
denote a spectral library with
3.
spectral signatures, each with I spectral bands, where it is
assumed that the number of spectral signatures is much larger
(m » I ) . Assuming that matrix 1It1 x n holds n observed spectral vectors
than the number of bands
y
. . . ,Y and is given by ==
[y 1 ,
E
n]
S
AS
+ N,
[SI , . . . , s n ] E lItm x n == [n l , . . . , n ] n
==
tion matrix and N
noise. To be physically meaningful are subj ect to
constraint
(1)
is the abundance frac
E
[4] ,
1It1 x n
is the additive
abundance fractions
nonnegativity constraint (ANC) and sum-to-one
(ASC). Thus, the abundance fractions estimation
problem can formulated as follow s :
s
where notation
S ;:::
In denote a
II (oj II F
x
m
and
1
x
n
1
shows the schematic of the proposed
In/321
of Algorithm l .
blocks of
block is
1 024,
32
x
I m/321
The grid contains
2 x
32 and the number of threads for each
thus all pixels in the hyperspectral image are
processed in parallel. The algorithm is constructed by a chain of kernels .
The first kernel computes the matrix
of Algorithm
1 ).
H
(line
2
In order to minimize the number of global
32
x
32,
which is the size of the block, and trans
ferred, step by step, to the shared memory. The result for each
(2)
stands for Frobenius norm, the in
° is to be understood componentwise,
1
9
and line
blocks of
=
subj ect to :
global memory. Fig .
implementation for the multiplication kernel used in line
memory accesses, matrices A and Y are partitioned into sub
�2 I I Y - M S I I } S ;::: 0 , I;z:,S 1;,
min
equality
SUNSAL it is highly parallelizable, since all calculations can be performed in a pixel-by-pixel basis. The first step is to map the hyperspectral image and the spectral library into GPU
y where
SUNSAL IMPLEMENTATION ON GPU
element of
H is the sum of several partial products. Inside the R (see line 8 of Algo
loop, another kernel compute matrix rithm
1 ).
This kernel launches as many threads as elements
R, where each thread computes an element H + p, (V k + D k ) and stores the result in the global
1m,
and
that are present in
column vectors filled of
1 's,
of the
respectively. D u e t o the fact that only a few signatures con
memory. The kernel that computes the abundance estimates,
S contains many zero values, which means that it is sparse.
on two operations: first it computes the product of matrices
tained in A are likely to contribute to the observed spectra Y,
on each iteration (line
and
9
of Algorithm
1 ),
is decomposed
R, followed by the addition of matrix cl � .
G
The scheme
herein used is similar to the first kernel used to compute ma
2.1. SUNSAL method The optimization problem
S,
trix
(2)
can be posed in the framework
of convex optimization, where the convex term promotes the
1447
H
whereas the kernels to update
V
same rule has the kernel used to compute updated by analyzing if each element of
and
R.
D
follow the
Finally,
V is negative.
V
is
Global Memory Matrix A
[m
x
L]
;
[L I I I I [ Matrix Y x
n
Table 1. Processing times
]
Cuprite dataset
I I I I
Matrix H m x n]
Shared Memory
342).
(l
=
1 88)
n
=
,,
,
I
I
I
1rn/ 3 21
:Block 2
:
:
2 -
--
:
----
:
1000
47750
n
=
1 0000
n
=
=
0 . 845
19.333
76.323
0.036
0 . 1 24
0.465
23.4
155.9
1 64. 1
,
': ::-1::::]:::::[ ::::::: : :[:::: - -- -1 [ 3 f � ----1- -- r r r
I-o!=f=+,�
(m
Sequential (CPU)
,,
,
seconds) and speedup for the
and for the USGS library
Parallel (GTX680) Speedup
,,
(l0 3
----
:
en
1) u
§ -g
- ---
::l
0 . 5 ,------�---�---, 0.4
�
"2 0. 3
O:i
E
.� 0.2
Fig. 1 . Illustration o f parallel SUNSAL in the GPU : multipli
1)
-5
cation kernel.
4-<
0 0. 1
