GENERALIZED BILINEAR MODEL BASED NONLINEAR UNMIXINGUSING SEMI-NONNEGATIVE MATRIX FACTORIZATION
Naoto Yokoya,1 Jocelyn Chanussot,2 and Akira Iwasaki3
1Department of Aeronautics and Astronautics, The University of Tokyo, Japan
GIPSA-lab, Signal and Image Dept., Grenoble Institute of Technology, Grenoble, France
3
Research Center for Advanced Science and Technology, The University of Tokyo, Japan

2ABSTRACT
Nonlinear spectral mixing models have recently been receiving attention in hyperspectral image processing. This work
presents a novel optimization method for nonlinear unmixing
based on a generalized bilinear model (GBM), which considers second-order scattering effects. Semi-nonnegative matrix factorization is used for optimization to process a whole
image in a matrix form. The proposed method is applied
to an airborne hyperspectral image with many endmembers
and shows good performance both in unmixing quality andcomputational cost with simple implementation. The effectof endmember extraction on nonlinear unmixing is investigated and the impact of the nonlinearity on abundance maps
is demonstrated.
Index Terms— hyperspectral image, nonlinear unmixing, generalized bilinear model, semi-nonnegative matrix factorization

with three endmembers [7, 8]. When applied to larger images in an unsupervised manner with many endmembers, theoptimization process becomes more challenging.
In this work, we present a novel optimization method for
the GBM and apply it to a real HS image with many endmembers. The effect of endmember extraction on nonlinear
unmixing is investigated and the impact of the second-orderscattering effects on abundance maps is visualized.
2. GBM UNMIXING VIA SEMI-NONNEGATIVE
MATRIX FACTORIZATION
2.1. Generalized bilinear model
The BMM takes account of second-order photon interactions
between different D endmembers as additional terms in theLMM assuming that third or higher order interactions are negligible [4]. In the BMM, the observed L-spectrum of a single
pixel y ∈ RL×1 is given by:

1. INTRODUCTION
Spectral unmixing is an important task for hyperspectral (HS)image interpretation. Many researchers have worked on thisproblem with a linear mixture model (LMM) that assumes
that an observed spectrum is a linear combination of several
endmember spectra. The LMM is a simplified spectral mixture model that considers only first-order scattered photons
by neglecting multiple photon interactions. Although LMM
based unmixing methods can retrieve physically meaningful
results, nonlinearity in spectral mixture model is pointed out
by some works [1]-[4]. In recent years, nonlinear unmixing
for HS images is receiving attention in remote sensing image
exploitation. Nonlinear spectral mixing occurs due to multiple reflection and transmission from surface [2]. The bilinear
mixture model (BMM), which considers second-order scattering of photons between two distinct materials, has beenstudied by several groups [5]-[7]. The generalized bilinear
model (GBM) introduces an effective means to deal with theunderlying assumptions in the BMM [7]. The GBM method
was applied to small images of synthetic and real HS datay = Ea +D−1∑

D∑bi,j ei ⊙ ej + n,(1)i=1 j=i+1

where E ∈ RL×D is the endmember matrix with the i th column vector ei ∈ RL×1 representing the i th endmember spectrum, a ∈ RD×1 is the abundance vector, bi,j is the interaction abundance between the i th and j th endmembers, ⊙ is theHadamard (element-wise product) operation, and n ∈ RL×1
is the additive noise. On the right side, the first term denotes
the linear mixing and the second term represents the secondorder interaction, which is assumed to be a linear combination of bilinear endmember spectra, i.e., the bilinear mixing.
From a physical perspective, the GBM introduces the nonlinear mixing coefficient ci,j as bi,j = ci,j ai aj and assumes thefollowing constraints:ai ≥ 0 ∀i ∈ 1, ..., D andD∑ai = 1,i=1

0 ≤ ci,j ≤ 1 ∀i ∈ 1, ..., D − 1 ∀j ∈ i + 1, ..., D.

(2)

When the endmembers are known, the GBM unmixing turns
to the optimization of the abundance (a) and the interaction
coefficient (c) under the constraints of (2). Several optimization methods are proposed in [8].
2.2. Semi-nonnegative matrix factorization for GBM
In this work, the new optimization method based on seminonnegative matrix factorization (Semi-NMF) is introduced
to speed up the optimization process of a whole image in amatrix form. The observed HS image can be reshaped as amatrix form Y ∈ RL×P with P representing the number ofpixels. The BMM for the whole image is given in a matrix
form by
Y = EA + MB + N(3)
where A ∈ RD×P is the abundance matrix, M ∈ RL×D(D−1)/2
is the bilinear endmember matrix, B ∈ RD(D−1)/2×P is theinteraction abundance matrix, and N ∈ RL×P is the noisematrix. The GBM unmixing becomes the following minimization with respect to A and B:
minimize ∥Y − EA − MB∥2F
subject to A ⪰ 0,D∑Ai,k = 1, 0 ⪯ B ⪯ A∗

(4)
(5)

i=1

where A∗(i,j),k = Ai,k Aj,k (∀k ∈ 1, ..., P ) and the operator∥ · ∥F denotes the Frobenius norm. By introducing Y1 =Y − MB and Y2 = Y − EA, (3) is written as follows
Y1 = EA + N,
Y2 = MB + N.

(6)
(7)

Separately minimizing ∥Y1 −EA∥2F and ∥Y2 −MB∥2F can be
an approximation of the solution of the original problem (4)performed by an alternating optimization algorithm. Owing
to physical constraints, all components of E, M, A, and B
are nonnegative. Therefore, the minimization of (4) can besolved by Semi-NMF that factorizes a non-restricted matrix
X into a non-restricted matrix F and a nonnegative matrix GTas X ≈ FGT [9]. Semi-NMF optimization is guaranteed toconverge to a local optimum with the alternative update rules.
With E given and M calculated from E, the GBM unmixing
is solved by the following update rules for A and B:√TTTTT+−−+AT ←AT .∗ ((YT(8)1 E) +A (E E) )./((Y1 E) +A (E E) )√TTTTT+−−+BT ←BT .∗ ((YT(9)2 M) +B (M M) )./((Y2 M) +B (M M) )where .∗ and ./ denote elementwise multiplication and division. (C)+ and (C)− are the positive and negative parts of amatrix C defined as C+ = (|C| + C)/2, C− = (|C| − C)/2.
The GBM unmixing via Semi-NMF is as follows.
Algorithm: Generalized bilinear model based nonlinear unmixing via Semi-NMF

Input: Hyperspectral data Y ∈ RL×P and endmember matrix E ∈ RL×D .
Output: Abundance matrix A ∈ RD×P and interaction abundance matrix B ∈ RD(D−1)/2×P .
Step 1: A is initialized by the fully constrained least square(FCLS) method [10] based on the LMM.
Step 2: A∗ is calculated and B is set as δ × A∗ with smallvalue of δ.
Step 3: A and B are alternately updated by (8) and (9). If any
element of B exceeds that of A∗ , it is replaced by that of A∗ .
To satisfy the abundance sum-to-one constraint, the methodfrom [10] is adopted.
3. EXPERIMENTThe GBM via Semi-NMF is applied to an airborne HS datataken over vegetation because multiple scattering is a common phenomenon in vegetation areas [1]-[4]. The dataset
was acquired by a Compact Airborne Spectrographic Imager
3 (CASI-3) taken over pasture area in Hokkaido, Japan, on
June 19, 2009. We selected a 150×150 pixel size image with68 spectral channels over 410-1070 nm and 1 m ground sampling distance (GSD). The data initially measured as radiancewere converted into reflectance.
Vertex component analysis (VCA) [11], which is one of
the most common convex-geometry-based endmember extraction methods with the pure pixel assumption, is used forendmember extraction. The performance of the proposed
method is numerically evaluated by the reconstruction error(RE), the spectral angle mapper (SAM), and the computational time comparing with FCLS and the gradient descentalgorithm (GDA), which showed good results for the GBM
unmixing applied to a real HS image in [8]. RE and SAM arecommonly used to evaluate the performance of the unmixing
procedure for real datasets. They are calculated using theobserved (yk ) and estimated (ŷk ) spectra as follows:vuP
u 1 ∑∥ŷk − yk ∥2 ,(10)
RE = tLPk=11 ∑arccos
PP

SAM =k=1(yk · ŷk∥yk ∥∥ŷk ∥).

(11)

By changing the number of endmembers extracted by VCA,
we examined the second-order interaction term and the residual errors of the GBM. The impact of these values at a single
pixel is calculated by root-sum-square (RSS). After the endmember extraction, the abundance maps and the impact ofbilinear mixing on these abundance maps are estimated.
4. RESULTS AND DISCUSSION
Table 1 shows the comparison of the RE and SAM values,
and the computational cost relative to different four numbers

Table 1. Comparison of RE (×10−3 ) and SAM (in degree) values, and computational costs (in sec).

FCLS
GDA
Semi-NMF

RE
10.562
10.54810.264

D=6
SAM
2.695
2.692
2.640

Time
20.5
122.2
55.4

RE
6.320
6.318
6.258

D = 11
SAM
1.824
1.823
1.816

of endmembers, i.e., D = 6, 11, 16, 21. The CPU used isIntel(R) Core(TM) i7 CPU 2.80 GHz, with a memory capacity of 16 GB. The Semi-NMF method shows the best performance in RE and SAM for many cases of the number ofendmembers. In addition, the computational time is smaller
than GDA. Especially, even when the number of endmembers
becomes large, the Semi-NMF method converges to good local optima with reasonable computational cost owing to thesimple update rules in a matrix form. This result indicatesthat Semi-NMF is a good candidate for the optimization ofthe GBM unmixing regarding its convergence and easy implementation.
Fig. 1 shows the RSS maps of the second-order interaction term and the unmixing residual errors. As the numberof endmembers increases (6, 11, 16, and 21), the two maps
converge. We concluded 16 is enough for the number of endmembers. In spite of the large number of endmembers, thebilinear mixing effect clearly appears in areas B and C in Fig.
1, i.e., the boundary between grass and tree, and the narrow
path. With D = 6, area A in Fig. 1 corresponding to grassshows the high bilinear effect but the residual error still remains high in this area. It significantly decreases after thegrass is detected as the endmember with D = 11. It means
the endmember extraction has larger influence on the residualerrors of unmixing than the bilinear mixing effect. Therefore,
to discuss the nonlinear mixing, the accurate endmember extraction is the major premise.
Since VCA is sensitive to the outlier, it detects some endmembers that are only present in very specific regions. To
visualize the abundance map and the bilinear mixing effect,
we manually eliminated such kind of specific endmembers
and re-processed the GBM unmixing with 10 endmembers
as shown in Fig. 2 (a). In Fig. 2 (b), the colored circles
indicate the locations of these endmembers, with each color
corresponding to that of spectrum in Fig. 2 (a). The endmembers are labeled by visual judgment considering the locationsof endmembers and comparing the abundance maps and theRGB color image. Fig. 2 (c) and (d) show the abundancemaps for these ten materials and their differences with theLMM, which demonstrate the impacts of the second-order interaction on the abundance maps. In some regions, the abundance fractions changed more than 2% from the LMM because of the multiple scattering effect.

Time
57.4
223.1
106.6

RE
5.025
5.048
4.975

D=6A

D = 16SAM
1.488
1.479
1.483

D=11Time
95.2
670.6
169.7

D = 21SAM
1.433
1.431
1.429

RE
4.738
4.721
4.701

D=16Time
143.5
785.7
250.4

D=210.12

B
C

00.2

0Fig. 1. RSS maps of second-order interaction terms (top) andunmixing residual errors (bottom).

5. CONCLUSION
This work presented a novel optimization method based on
Semi-NMF for the GBM unmixing and demonstrated its application to an airborne real HS image taken over a vegetation scene with many endmembers. Semi-NMF enables theoptimization of the GBM unmixing for processing a whole
image in a matrix form, with simple update rules. The proposed method showed useful results both in unmixing quality
and computational cost. By examining the effect of endmember extraction on nonlinear unmixing, the residual errors ofunmixing are more influenced by the endmember extraction
than the bilinear spectral mixing effect. The impact of thesecond-order scattering effects on abundance maps is visualized. With accurate endmember extraction, Semi-NMF based
GBM has a possibility to deal with spectrally complex HS
images, thereby contributing to the practical use of nonlinear
unmixing.

6. ACKNOWLEDGMENT
The authors would like to thank Prof. J. Tourneret for providing the GDA code. This project is partly carried out under thecontracted with AIST (National Institute of Advanced Industrial Science and Technology).

(a)(c) tree

0.7
0.6
0.5
Reflectancetree (shadow)

grass 1

grass 2

grass 3treetree (shadow)
grass 1
grass 2
grass 3
grass 4
soil 1soil 2road 1road 20.40.31

grass 4

soil 1soil 2road 1road 2

0.2

00.1
0400

500

600

700
800
900
Wavelength (nm)1000

1100(b)(d) treetree (shadow)

grass 1

grass 2

grass 3-2×102

grass 4

soil 1soil 2road 1road 2-2

Fig. 2. (a) Endmember spectra, (b) location of endmembers, (c) abundance maps, and (d) differences of abundance maps withLMM.

7. REFERENCES[1] D.A. Roberts, M.O. Smith, and J.B. Adams, “Greenvegetation, nonphotosynthetic vegetation, and soils in
AVIRIS data,” Remote Sens. Environ., vol. 44, pp. 117–126, 1993.
[2] C. C. Borel and S. A. W. Gerstl, “Nonlinear spectral
mixing models for vegetative and soil surfaces,” RemoteSens. Environ., vol. 47, pp. 403–416, 1994.
[3] L. Zhang, D. Li, Q. Tong, and L. Zheng, “Study of thespectral mixture model of soil and vegetation in PoYang
Lake area, China,” Int. J. Remote Sens., vol. 19, no. 11,
pp. 2077–2084, 1998.
[4] B. Somers, K. Cools, S. Delalieux, J. Stuckens, D. V. derZande, W. W. Verstraeten, and P. Coppin, “Nonlinear
hyperspectral mixture analysis for tree cover estimatesin orchards,” Remote Sens. Environ., vol. 113, no. 6, pp.
1183–1193, Jun. 2009.
[5] W. Fan, B. Hu, J. Miller, and M. Li, “Comparative
study between a new nonlinear model and common linear model for analysing laboratory simulated-forest hyperspectral data,” Int. J. Remote Sens., vol. 30, no. 11,
pp. 2951–2962, Jun. 2009.

[6] J. M. Nascimento and J. M. Bioucas-Dias, “Nonlinear
mixture model for hyperspectral unmixing,” in Proc.
SPIE, 2009, vol. 7477, p. 74770I.
[7] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret,“Nonlinear unmixing of hyperspectral images using ageneralized bilinear model,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 11, pp. 4153–4162, Nov. 2011.
[8] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret,“Unmixing hyperspectral images using the generalizedbilinear model,” IEEE IGARSS, pp. 1886–1889, 2011.
[9] C. Ding, T. Li, and M. I. Jordan, “Convex and seminonnegative matrix factorization,” IEEE Trans. Patten
Anal. Mach. Intell., vol. 32, no. 1, pp. 45–55, 2010.
[10] D. C. Heinz and C.-I. Chang, “Fully constrained leastsquares linear spectral mixture analysis method for material quantification in hyperspectral imagery,” IEEE
Trans. Geosci. Remote Sens., vol. 39, no. 3, pp. 529–545, Mar. 2001.
[11] J. M. Nascimento and J. M. Bioucas-Dias, “Vertex component analysis: A fast algorithm to unmix hyperspectral data,” IEEE Trans. Geosci. Remote Sens., vol. 43,
no. 4, pp. 898–910, Apr. 2005.