<figure><img src="image_1.png" /></figure>

HAL
open science

# Using Low-Rank Representation of Abundance Maps and Nonnegative Tensor Factorization for Hyperspectral Nonlinear Unmixing

Lianru Gao, Zhicheng Wang, Lina Zhuang, Haoyang Yu, Bing Zhang, Jocelyn Chanussot

<figure><img src="image_2.png" /></figure>

## To cite this version:

Lianru Gao, Zhicheng Wang, Lina Zhuang, Haoyang Yu, Bing Zhang, et al. Using Low-Rank Representation of Abundance Maps and Nonnegative Tensor Factorization for Hyperspectral Non-linear Unmixing. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60, pp.5504017. 10.1109/TGRS.2021.3065990. hal-03932934

HAL Id: hal-03932934

https://hal.science/hal-03932934v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Using Low-rank Representation of Abundance Maps and Nonnegative Tensor Factorization for Hyperspectral Nonlinear Unmixing

Lianru Gao, Senior Member, IEEE, Zhicheng Wang, Lina Zhuang, Member, IEEE, Haoyang Yu, Member, IEEE, Bing Zhang, Fellow, IEEE, and Jocelyn Chanussot, Fellow, IEEE

**Abstract**—Tensor-based methods have been widely studied to attack inverse problems in hyperspectral imaging, since a hyperspectral image cube can be naturally represented as a third-order tensor, which can perfectly retain the spatial information in the image. In this paper, we extend the linear tensor method to the nonlinear tensor method, and propose a nonlinear low-rank tensor unmixing algorithm to solve the generalized bilinear model (GBM). Specifically, the linear and nonlinear parts of the GBM can both be expressed as tensors. Furthermore, the low-rank structures of abundance maps and nonlinear interaction abundance maps are exploited by minimizing their nuclear norm, thus taking full advantage of the high spatial correlation in hyperspectral images. Synthetic and real-data experiments show that the low rank of abundance maps and nonlinear interaction abundance maps exploited in our method can improve the performance of the nonlinear unmixing. **A MATLAB demo of this work will be available at https://github.com/LinaZhuang for the sake of reproducibility.**

**Index Terms**—Hyperspectral Image, Tensor Decomposition, Low Rank, Nonlinear Unmixing.

## I. INTRODUCTION

HYPERSPECTRAL remote sensing imaging technology has developed rapidly in the past 20 years and is playing a key role in Earth observation. Hyperspectral images (HSIs) are acquired by hyperspectral cameras (HSCs) that have hundreds or thousands of spectral bands and a rich spectral resolution ranging from 0.4 ~ 2.5µm [1]–[3]. Because of the low spatial resolution of HSCs, microscopic material mixing, and the multiple scattering that occurs within scenes [2], the spectral signals in HSIs consist of mixtures of spectra of

different materials. However, one benefit of the high spectral resolution of HSIs is that the problem of mixed pixels can potentially be solved. Hyperspectral unmixing (HU) aims to obtain the basic components of the image, called endmembers, and their corresponding proportions, called abundances [1]. To solve this problem, there are two main approaches that can be used. The linear mixing model (LMM) assumes that photons interact with only one material before reaching the sensor. This linear model has been widely developed and tested in recent years [3]–[7], but is only suitable for relatively simple scenarios. The second main method, hyperspectral nonlinear unmixing (HNU), consists of various nonlinear mixing models (NLMMs) that consider the multiple or infinite reflections of photons. This method has been shown to be more appropriate for application to real scenes. This article focuses on HNU.

The NLMMs can be divided into two main categories: bilinear mixing models (BMMs) and high-order mixing models. The former group includes, for example, the Fan model [8], the generalized bilinear model (GBM) [9], and the linear-quadratic mixing model (LQM) [10]. The latter group includes the multilinear mixing model (MLM) [11], the *p*-linear model [12], multiharmonic postnonlinear mixing model (MHPNMM) [13], etc. The GBM is one of the most popular of the BMMs since it is more suitable than the Fan model (FM) for modeling scenes, and it is also a generalization of both the LMM and FM [3]. Meanwhile, the GBM can more flexibly quantify the strength of the contributions of different bilinear components than the polynomial post nonlinear model (PPNM) [3], [14]. In this work, we address the spectral unmixing problem based on the GBM.

The Bayesian and gradient descent algorithm (GDA) are widely used to solve optimization problems. They have been used to solve the GBM in [15]: in this case the image was unmixed pixel by pixel and similar results were obtained. To deal with the high computational cost of pixel-based unmixing, the semi-nonnegative matrix factorization (semi-NMF) algorithm has been proposed for accelerating the optimization of whole images in matrix form [16]. A bound projected optimal gradient method that transforms the GBM into a least-squares problem was also proposed in [17]. In addition, a multi-task learning (MTL) method that assumes that the linear part of the GBM is the task to be solved and the nonlinear part is a secondary task was proposed in [18]. However, all of these methods are based on the hyperspectral image matrix, which means that the original spatial structure is destroyed

This work was supported in part by the National Natural Science Foundation of China under Grant 42030111, Grant 41722108, and Grant 42001287 and in part by the AXA Research Fund. (Corresponding author: Lina Zhuang.)

L. Gao is with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China (e-mail: gaolr@aircas.ac.cn).

Z. Wang and B. Zhang are with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China, and also with the University of Chinese Academy of Sciences, Beijing 100049, China (e-mail: wangzc@radi.ac.cn; zb@radi.ac.cn).

L. Zhuang is with the Department of Mathematics, Hong Kong Baptist University, Hong Kong (e-mail: linazhuang@qq.com).

H. Yu is with the Center of Hyperspectral Imaging in Remote Sensing (CHIRS), Information Science and Technology College, Dalian Maritime University, Dalian, 116026, China (e-mail: yuhy@dlmu.edu.cn).

Jocelyn Chanussot is with the Université Grenoble Alpes, INRIA, CNRS, Grenoble INP, LJK, 38000 Grenoble, France, and also with Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China (e-mail: jocelyn@hi.is).
and spatial information is lost. As the HSI can be naturally expressed as a third-order tensor, which combines spectral and spatial information naturally, tensor decomposition-based methods have been applied to hyperspectral image analysis using canonical polyadic (CP) decomposition [19], Tucker decomposition [20], and other methods.

The low rank of abundances has been exploited to solve the spectral unmixing problem [21], [22]. Abundances are organized as a two-dimensional matrix $\mathbf{A} = [\mathbf{a}_1, \dots, \mathbf{a}_N] \in \mathbb{R}^{R \times N}$, whose columns contain $N$ abundance vectors corresponding to $N$ pixels and rows which are abundance maps corresponding to $R$ endmembers (or atoms in a dictionary). In order to take advantage of the correlation among abundances, the low-rank structure of abundance matrix $\mathbf{A}$ was exploited in [14], [18], [23]–[29], by minimizing the nuclear norm of matrix $\mathbf{A}$ [23], [29] in optimization problems.

Apart from investigating the low-rank abundance matrix, the low-rank structure of abundance maps was first discussed in [30] in relation to solving the linear-mixing problem. Abundances corresponding to a single endmember were organized as an abundance map, which is a two-dimensional matrix of size $n_{row} \times n_{col}$ (where $n_{row}$ and $n_{col}$ are the number of rows and columns in the HSI). Due to the high spatial correlation in the HSIs and the sparse distribution of endmembers, abundance maps can be approximated using low-rank representations. In [30] the low rank of abundance maps was enforced by matrix factorization: i.e., the abundance map was represented as matrix that was the product of two smaller matrices and whose size determined the rank of the abundance map. Instead of using matrix factorization, we intend to enforce the low rank of abundance maps by minimizing the nuclear norm of abundance maps, which does not require prior knowledge about the rank of the abundance maps.

This paper addresses the nonlinear unmixing problem by taking advantage of the low rank of the abundance maps and nonlinear interaction maps. We introduce a new nonlinear unmixing method based on a tensor decomposition algorithm and low-rank representation to solve the GBM. The main contributions of this work can be summarized as follows.

* We proposed a new unmixing method based on the GBM model. Instead of using a matrix format, GBM is rewritten in a tensor format, which is a good representation for embedding the inherent spectral-spatial structure of HSI. The tensor format allows us to address the nonlinear unmixing problem as a nonnegative tensor factorization problem and to exploit the spatial correlation of abundance maps.

* The idea of the low-rank representation of abundance maps is extended to the nonlinear mixing components - i.e., the low-rank representation of nonlinear interaction maps. The spatial pattern of nonlinear mixing components is discussed in detail in Section II-C-2.

The rest of the paper is organized as follows. Section II introduces the GBM model and the proposed tensor-based algorithm. Section III and Section IV report results for the synthetic data and real dataset, respectively. Section V concludes the paper.

## II. PROBLEM FORMULATION AND METHOD

In this section, we first introduce the generalized bilinear model (GBM). The proposed nonnegative tensor factorization-based method for estimating abundances is then described.

### A. Notation and Definitions

In this subsection, we introduce the notation and definitions used in the paper. An nth-order tensor is identified using Euler-cript letters - e.g., $Q \in \mathbb{R}^{I_1 \times I_2 \times \dots \times I_n \times \dots \times I_N}$, with the $I_n$ is the size of the corresponding dimension $n$. Hence, an HSI can be naturally represented as a third-order tensor, $T \in \mathbb{R}^{I_1 \times I_2 \times I_3}$, which consists of $I_1 \times I_2$ pixels and $I_3$ spectral bands. Three further definitions related to tensors are given as follows.

**Definition 1:** The dimension of a tensor is called the mode: $Q \in \mathbb{R}^{I_1 \times I_2 \times \dots \times I_n \times \dots \times I_N}$ has $N$ modes. For a third-order tensor $T \in \mathbb{R}^{I_1 \times I_2 \times I_3}$, by fixing one mode, we can obtain the corresponding sub-arrays, called slices - e.g., $T_{i,:;n}$.

**Definition 2:** The mode-n product of a tensor $Q \in \mathbb{R}^{I_1 \times I_2 \times \dots \times I_N}$ by a matrix $\mathbf{X} \in \mathbb{R}^{J_n \times I_n}$ is a tensor $\mathcal{G} \in \mathbb{R}^{I_1 \times \dots \times I_{n-1} \times J_n \times I_{n+1} \times \dots \times I_N}$, denoted as

$$ \mathcal{G} = Q \times_n \mathbf{X}, $$

where each entry of $\mathcal{G}$ is defined as the sum of products of corresponding entries in $Q$ and $\mathbf{X}$:

$$ \mathcal{G}(i_1, \dots, i_{n-1}, j_n, i_{n+1}, \dots, i_N) = \sum_{i_n} \mathcal{Q}(i_1, \dots, i_N) \cdot \mathbf{X}(j_n, i_n). $$

**Definition 3:** Given a matrix $\mathbf{A} \in \mathbb{R}^{k_1 \times k_2}$ and vector $\mathbf{c} \in \mathbb{R}^{l_1}$, their outer product, denoted as $\mathbf{A} \circ \mathbf{c}$, is a tensor with dimensions $(k_1, k_2, l_1)$ and entries $(\mathbf{A} \circ \mathbf{c})_{i_1, i_2, j_1} = \mathbf{A}_{i_1, i_2} \mathbf{c}_{j_1}$.

### B. Spectral Nonlinear Model: GBM

Bilinear mixture models (BMMs) are based on considerations of the second-order interactions between different endmembers. These models can overcome the inherent limits of the linear model and can extract complex information from the scene to improve the unmixing results. By considering bilinear interactions as additional endmembers, a pixel $\mathbf{y} \in \mathbb{R}^{L \times 1}$ with $L$ spectral bands can be expressed as follows:

$$ \mathbf{y} = \mathbf{C}\mathbf{a} + \sum_{i=1}^{R-1} \sum_{j=i+1}^{R} b_{i,j}\mathbf{c}_i \odot \mathbf{c}_j + \mathbf{n}, \quad (1) $$

where $\mathbf{C} = [\mathbf{c}_1, \mathbf{c}_2, \dots, \mathbf{c}_R] \in \mathbb{R}^{L \times R}$, $\mathbf{a} = [a_1, a_2, \dots, a_R]^T \in \mathbb{R}^{R \times 1}$, and $\mathbf{n} \in \mathbb{R}^{L \times 1}$ represent the mixing matrix containing the spectral signatures of $R$ endmembers, the fractional abundance vector, and the white Gaussian noise, respectively. The nonlinear coefficient $b_{i,j}$ controls the nonlinear interaction between the materials, and $\odot$ is a Hadamard product operation.

To satisfy the physical assumptions and overcome the limitations of the FM [3], the GBM redefines the parameter $b_{i,j}$ as $b_{i,j} = \gamma_{i,j}a_ia_j$. Meanwhile, the abundance non-negativity
constraint (ANC) and the abundance sum-to-one constraint (ASC) are satisfied as follows:

$$ 
\begin{aligned}
& a_i \ge 0, \sum_{i=1}^{R} a_i = 1, \\
& 0 < \gamma_{i,j} < 1, \forall i < j, \\
& \gamma_{i,j} = 0, \forall i \ge j.
\end{aligned} 
\quad (2) 
$$

The spectral mixing model for N pixels can be written in matrix as:

$$ \mathbf{Y} = \mathbf{C}\mathbf{A} + \mathbf{M}\mathbf{B} + \mathbf{N}, \quad (3) $$

where $\mathbf{Y} = [\mathbf{y}_1, \mathbf{y}_2, \dots, \mathbf{y}_N] \in \mathbb{R}^{L \times N}$, $\mathbf{A} = [\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_N] \in \mathbb{R}^{R \times N}$, $\mathbf{M} \in \mathbb{R}^{L \times R(R-1)/2}$, $\mathbf{B} \in \mathbb{R}^{R(R-1)/2 \times N}$, and $\mathbf{N} \in \mathbb{R}^{L \times N}$ represent the hyperspectral image matrix, the fractional abundance matrix with N abundance vectors (the columns of A), the bilinear interaction endmember matrix, the nonlinear interaction abundance matrix, and the white Gaussian noise matrix, respectively.

This work aims to solve a model-based supervised unmixing problem: that is to estimate the abundances, A, and nonlinear coefficients, B, given the spectral signatures of the endmembers, C, which are known as a prior.

### C. Nonlinear Unmixing Based on Nonnegative Tensor Factorization

1) Tensor Framework of the GBM: Traditional nonlinear unmixing methods, such as the GDA and semi-NMF, transform the hyperspectral image cube into a two-dimensional matrix for processing, thus destroying the internal spatial structure of the data and resulting in poor abundance inversion. However, given that the hyperspectral images can be naturally represented as a third-order tensor, we propose a hyperspectral nonlinear unmixing algorithm based on tensor representation for the original hyperspectral image cube. The hyperspectral image cube $\mathcal{Y} \in \mathbb{R}^{n_{row} \times n_{col} \times L}$ can be expressed in the following format:

$$ \mathcal{Y} = \mathcal{A} \times_3 \mathbf{C} + \mathcal{B} \times_3 \mathbf{M} + \mathcal{N}, \quad (4) $$

where $\mathcal{A} \in \mathbb{R}^{n_{row} \times n_{col} \times R}$, $\mathcal{B} \in \mathbb{R}^{n_{row} \times n_{col} \times R(R-1)/2}$, and $\mathcal{N} \in \mathbb{R}^{n_{row} \times n_{col} \times L}$ denote the abundance cube containing R endmembers, the nonlinear interaction abundance cube, and the white Gaussian noise cube, respectively.

The tensor model (4) is equivalent to the matrix model (3), however, using a tensor form, enables us to remain the spatial structure of abundances in the model, so that its spatial structure can be exploited by adding regularizations.

To better represent the structure of abundance maps, mixing model (4) can be equivalently written as

$$ \mathcal{Y} = \sum_{i=1}^{R} \mathcal{A}_{:,i} \circ \mathbf{c}_i + \sum_{j=1}^{R(R-1)/2} \mathcal{B}_{:,j} \circ \mathbf{m}_j + \mathcal{N}, \quad (5) $$

where $\mathcal{A}_{:,i} \in \mathbb{R}^{n_{row} \times n_{col}}$, $\mathbf{c}_i \in \mathbb{R}^{L \times 1}$, $\mathcal{B}_{:,j} \in \mathbb{R}^{n_{row} \times n_{col}}$, and $\mathbf{m}_j \in \mathbb{R}^{L \times 1}$ denote the $i$th abundance slice, the $j$th endmember vector, the $k$th interaction abundance slice, and the $l$th interaction endmember vector, respectively. Model (5) is depicted in Fig. 1.

2) Low-rankness of Interaction Abundance Maps: As geographic data, the materials in HSIs tend to be spatially dependent. Spatial dependence means that things that are spatially close together tend to be more closely related than things that are far apart. Thus, the spatial distribution of a single material/endmember tends to be aggregated instead of being purely random. This spatial distribution of an endmember enables us to approximate its abundance map with a low-rank matrix. The abundance map of a material was represented as a matrix that was the product of two smaller matrices and whose size determined the rank of the abundance map, which was firstly introduced for the linear unmixing model in [30]. In this paper, we further study the low-rank representation of interaction abundance maps. We exploit this spatial pattern of endmembers using a low-rank approximation for the interaction abundance map.

The low-rank structure of abundance maps and interaction abundance maps can be illustrated using an example. The nonlinear unmixing of a Cuprite image of size 250 × 191 pixels was carried out. Fig. 2 shows the abundance maps and interaction abundance maps estimated by the GDA [15] followed by an endmember estimation step (which uses vertex component analysis (VCA) [31]). To see the matrix structure of abundance maps, we performed singular value decomposition (SVD) of each abundance map $\mathcal{A}_{:,i}$ ($i = 1, \dots, R$) and each interaction abundance map $\mathcal{B}_{:,j}$ ($j = 1, \dots, R(R-1)/2$). Singular values are plotted in the final column of Fig. 2, where blue curves represent the cumulative probability distributions of singular values. When the cumulative probability curve reaches 95%, the corresponding dimension is marked using a pink dashed line. Take the first row of Fig. 2 as an abundance maps' example. It can be seen that 95% of the data variability can be represented well in a subspace with 126 dimensions. A low-rank representation of the abundance map of the endmember 'Sphene' is given in Fig. 2(b), and is a very good approximation of the original abundance map (as implied by the difference map in Fig. 2(c)). Furthermore, it can be seen that for interaction abundance maps shown in Fig. 2(m) and Fig. 2(q), the cumulative probability curves reaches 95% when the corresponding dimensions are 85 and 21, respectively, which are much smaller than the subspace dimensions in the abundance maps. The corresponding low-rank representation of the interaction abundance maps are shown Fig. 2(n) and Fig. 2(r), which indicate that the interaction abundance maps can be approximated well by low-rank matrices. We remark that the ranks of abundance maps is usually not as low as the rank of image matrix, $\mathbf{Y}$, which usually can be approximated using a subspace with dimension lower than 20. The rank of a abundance map depends on the spatial distribution of the corresponding endmember. As long as a abundance map is not a full-rank matrix, its low-rankness can be imposed in the objective function to obtain a better estimate of the abundance map. The lower the rank is, the better result of the estimated abundance map we can expect. It has been demonstrated in the work [30] that the use of low-rankness of abundance maps can improve linear unmixing performance. In Fig. 2, we justify that the ranks of interaction abundance maps are lower than the ranks of abundance maps, implying the use of low-rankness
<figure><img src="image_4.png" /><figcaption>Fig. 1. The generalized bilinear model using the tensor-based framework</figcaption></figure>

of interaction abundance maps can also improve nonlinear unmixing performance.

To take full advantage of the low-rank structure of abundance maps and interaction abundance maps, we propose a new nonlinear unmixing method based both on the **Low-Rank** representation of abundance maps and interaction abundance maps via **Nonnegative Tensor Factorization**, termed **LR-NTF**, which aims to solve the following optimization problem:

$$
\begin{align}
\underset{A_{:,;,i} \ge 0, A_{:,;,j}^* \ge B_{:,;,j} \ge 0}{\operatorname{arg min}} & \quad \frac{1}{2} \left\| \mathcal{Y} - \sum_{i=1}^{R} A_{:,;,i} \circ \mathbf{c}_i - \sum_{j=1}^{R(R-1)/2} B_{:,;,j} \circ \mathbf{m}_j \right\|^2 \nonumber \\
& \quad + \lambda_1 \left( \sum_{i=1}^{R} \| A_{:,;,i} \|_* \right) + \lambda_2 \left( \sum_{j=1}^{R(R-1)/2} \| B_{:,;,j} \|_* \right) \nonumber \\
\text{s.t.} & \quad \sum_{i=1}^{R} A_{:,;,i} = \mathbf{1}_{n_{row}} \mathbf{1}_{n_{col}}^T, \tag{6}
\end{align}
$$

where $A_{k_1,k_2,(i,j)}^* = A_{k_1,k_2,i} A_{k_1,k_2,j} (k_1 \in \{1, \dots, n_{row}\}, k_2 \in \{1, \dots, n_{col}\},$ and $(i,j) \in \{(i,j) | i < j, i \in \{1, \dots, R\}, j \in \{1, \dots, R\}\}$), and $\|\mathcal{X}\|_F^2$ denotes the Frobenius norm which returns the square root of the sum of the absolute squares of its elements. The symbol $\|\cdot\|_*$ represents the nuclear norm, and $\mathbf{1}_d$ represents a vector whose components are all one and whose dimension is given by its subscript. Let $\mathbf{P} = \mathbf{1}_{n_{row}} \mathbf{1}_{n_{col}}^T$. Abundance maps $A_{:,;,i}$ and $B_{:,;,j}$ are enforced to be low-rank by minimizing their nuclear norms. $\lambda_1 > 0$ and $\lambda_2 > 0$ are parameters of regularizations.

The optimization problem in (6) can be solved by optimization using alternating direction method of multipliers (ADMM) [32]. To use the ADMM, first (6) is converted into an equivalent form by introducing multiple auxiliary variables $\mathbf{V}_i$, $\mathbf{E}_j$ to replace $A_{:,;,i}$ ($i = 1, \dots, R$), $B_{:,;,j}$ ($j = 1, \dots, R(R-1)/2$). The formulation is as follows:

$$
\begin{align}
& \underset{A_{:,;,i} \ge 0, A_{:,;,j}^* \ge B_{:,;,j} \ge 0}{\operatorname{arg min}} \quad \frac{1}{2} \left\| \mathcal{Y} - \sum_{i=1}^{R} A_{:,;,i} \circ \mathbf{c}_i - \sum_{j=1}^{R(R-1)/2} B_{:,;,j} \circ \mathbf{m}_j \right\|_F^2 \nonumber \\
& \qquad + \lambda_1 \left( \sum_{i=1}^{R} \| \mathbf{V}_i \|_* \right) + \lambda_2 \left( \sum_{j=1}^{R(R-1)/2} \| \mathbf{E}_j \|_* \right) \nonumber \\
& \qquad \qquad \left\{ \begin{array}{@{}l@{}}
A_{:,;,i} = \mathbf{V}_i, i = 1, 2, \dots, R \\
B_{:,;,j} = \mathbf{E}_j, j = 1, 2, \dots, R(R-1)/2 \\
\sum_{i=1}^{R} A_{:,;,i} = \mathbf{1}_{n_{row}} \mathbf{1}_{n_{col}}^T
\end{array} \right. . \tag{7}
\end{align}
$$

By using the Lagrangian function, (7) can be reformulated

$$
\begin{align}
& \mathcal{L}(A_{:,;,i}, B_{:,;,j}, \mathbf{V}_i, \mathbf{E}_j, \mathbf{D}_i, \mathbf{H}_j, \mathbf{G}) = \nonumber \\
& \frac{1}{2} \left\| \mathcal{Y} - \sum_{i=1}^{R} A_{:,;,i} \circ \mathbf{c}_i - \sum_{j=1}^{R(R-1)/2} B_{:,;,j} \circ \mathbf{m}_j \right\|_F^2 \nonumber \\
& + \lambda_1 \left( \sum_{i=1}^{R} \| \mathbf{V}_i \|_* \right) \nonumber \\
& + \lambda_2 \left( \sum_{j=1}^{R(R-1)/2} \| \mathbf{E}_j \|_* \right) \nonumber \\
& + \frac{\mu}{2} \left( \sum_{i=1}^{R} \| \mathbf{V}_i - \mathbf{D}_i \|_F^2 \right) \nonumber \\
& + \frac{\mu}{2} \left( \sum_{j=1}^{R(R-1)/2} \| B_{:,;,j} - \mathbf{E}_j - \mathbf{H}_j \|_F^2 \right) \nonumber \\
& + \frac{\mu}{2} \left\| \sum_{i=1}^{R} A_{:,;,i} - \mathbf{P} - \mathbf{G} \right\|_F^2, \tag{8}
\end{align}
$$

where $\mathbf{D}_i$, $\mathbf{H}_j$ and $\mathbf{G}$ are Lagrange multipliers, and $\mu$ is the penalty parameter. The variables $A_{:,;,i}$, $B_{:,;,j}$, $\mathbf{V}_i$, $\mathbf{E}_j$, $\mathbf{D}_i$, $\mathbf{H}_j$, $\mathbf{G}$ were updated sequentially: this step is shown in **Algorithm 1**. The optimization details of the loss function are as given in the following:

Update $\mathcal{A}$. The optimization problem for $A_{:,;,i}$ is

$$
\begin{align}
A_{:,;,i}^{k+1} &= \arg \min_{A_{:,;,i}^k} \frac{1}{2} \left\| \mathcal{Y} - \sum_{i=1}^{R} A_{:,;,i}^k \circ \mathbf{c}_i - \sum_{j=1}^{R(R-1)/2} B_{:,;,j}^k \circ \mathbf{m}_j \right\|_F^2 \nonumber \\
&\quad + \frac{\mu}{2} \left\| A_{:,;,i}^k - \mathbf{V}_i^k - \mathbf{D}_i^k \right\|_F^2 + \frac{\mu}{2} \left\| \sum_{i=1}^{R} A_{:,;,i}^k - \mathbf{P} - \mathbf{G}^k \right\|_F^2 \nonumber \\
&= \arg \min_{A_{:,;,i}^k} \frac{1}{2} \sum_{b=1}^{L} \| \mathcal{O}_{:,;,b} - A_{:,;,i}^k c_{b_i} \|_F^2 + \frac{\mu}{2} \left( \sum_{i=1}^{R} \| A_{:,;,i}^k - \mathbf{V}_i^k - \mathbf{D}_i^k \|_F^2 \right) \nonumber \\
&\quad + \frac{\mu}{2} \left\| A_{:,;,i}^k + \tilde{\mathbf{A}} - \mathbf{P} - \mathbf{G}^k \right\|_F^2, \tag{9}
\end{align}
$$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Low-rank approximation of abundance maps and nonlinear interaction abundance maps in the Cuprite image</figcaption></figure>
where $\mathcal{O} = \mathcal{Y} - \sum_{i=1, \neg i}^{R} \mathcal{A}_{:,;,i}^{k} \circ \mathbf{c}_i - \sum_{j=1}^{R(R-1)/2} \mathcal{B}_{:,;,j}^{k} \circ \mathbf{m}_j \in \mathbb{R}^{n_{row} \times n_{col} \times L}$, and $\mathcal{O}_{:,;,b}$ is the $b$th slice. Meanwhile, $\tilde{\mathbf{A}} = \sum_{i=1, \neg i}^{R} \mathcal{A}_{:,;,i}^{k} \in \mathbb{R}^{n_{row} \times n_{col}}$ represents the sum of all abundance

maps apart from the $i$th slice, $\sum_{i=1, \neg i}^{R} a_i = a_1 + \cdots + a_{i-1} + a_{i+1} + \cdots + a_R$, and $\mathbf{c}_i = [c_{i_1}, c_{i_2}, \dots, c_{i_b}, \dots, c_{i_L}]^T \in \mathbb{R}^{L \times 1}$ is the $i$th endmember. Hence the solution for $\mathcal{A}_{:,;,i}$ can be derived as follows:

$$ \mathcal{A}_{:,;,i}^{k+1} = (\sum_{b=1}^{L} c_{b_i} c_{b_i}^T + 2\mu \mathbf{I})^{-1} (\sum_{b=1}^{B} \mathcal{O}_{:,;,b} c_{b_i}^T + \mu (\mathbf{V}_i^k + \mathbf{D}_i^k + \mathbf{P} + \mathbf{G}^k - \tilde{\mathbf{A}})). \quad (10) $$

Update $\mathcal{B}$. The optimization problem for $\mathcal{B}_{:,;,j}$ is

$$ \begin{aligned} \mathcal{B}_{:,;,j}^{k+1} &= \arg \min_{\mathcal{B}_{:,;,j}^{k}} \frac{1}{2} \left\| \mathcal{Y} - \sum_{i=1}^{R} \mathcal{A}_{:,;,i}^{k+1} \circ \mathbf{c}_i - \sum_{j=1}^{R(R-1)/2} \mathcal{B}_{:,;,j}^{k} \circ \mathbf{m}_j \right\|_F^2 \\ &\quad + \frac{\mu}{2} \| \mathcal{B}_{:,;,j}^{k} - \mathbf{E}_j^k - \mathbf{H}_j^k \|_F^2 \\ &= \frac{1}{2} \sum_{b=1}^{L} \| \mathcal{K}_{:,;,b} - \mathcal{B}_{:,;,j}^{k} m_{b_j} \|_F^2 + \frac{\mu}{2} \| \mathcal{B}_{:,;,j}^{k} - \mathbf{E}_j^k - \mathbf{H}_j^k \|_F^2, \end{aligned} \quad (11) $$

where $\mathcal{K} = \mathcal{Y} - \sum_{i=1}^{R} \mathcal{A}_{:,;,i}^{k} \circ \mathbf{c}_i - \sum_{j=1, \neg j}^{R(R-1)/2} \mathcal{B}_{:,;,j}^{k} \circ \mathbf{m}_j \in \mathbb{R}^{n_{row} \times n_{col} \times L}$, and $\mathcal{K}_{:,;,b}$ is the $b$th slice. Meanwhile, $\mathbf{m}_j = [m_{j_1}, m_{j_2}, \dots, m_{j_b}, \dots, m_{j_L}]^T \in \mathbb{R}^{L \times 1}$ is the $j$th interaction endmember. Hence the solution for $\mathcal{B}_{:,;,j}$ can be derived as follows:

$$ \mathcal{B}_{:,;,j}^{k+1} = (\sum_{b=1}^{L} m_{b_j} m_{b_j}^T + \mu \mathbf{I})^{-1} (\sum_{b=1}^{L} \mathcal{K}_{:,;,b} m_{b_j}^T + \mu (\mathbf{E}_j^k + \mathbf{H}_j^k)). \quad (12) $$

Update $\mathbf{V}$. The optimization problem for $\mathbf{V}_i$ is

$$ \begin{aligned} \mathbf{V}_i^{k+1} &= \arg \min_{\mathbf{V}_i^k} \lambda_1 \|\mathbf{V}_i^k\|_* + \frac{\mu}{2} \|\mathcal{A}_{:,;,i}^{k+1} - \mathbf{V}_i^k - \mathbf{D}_i^k\|_F^2 \\ &= \frac{1}{2} \|\tilde{\mathbf{V}}_i - \mathbf{V}_i^k\|_F^2 + \frac{\lambda_1}{\mu} \|\mathbf{V}_i^k\|_*; \end{aligned} \quad (13) $$

where $\tilde{\mathbf{V}}_i = \mathcal{A}_{:,;,i}^{k+1} - \mathbf{D}_i^k \in \mathbb{R}^{n_{row} \times n_{col}}$. Sub-problem (13) can be solved using singular value composition (SVD) of $\tilde{\mathbf{V}}_i$:

$$ [\mathbf{U}, \mathbf{S}, \mathbf{Z}] = \text{SVD}(\tilde{\mathbf{V}}_i), \quad (14) $$

where $\mathbf{U}$ and $\mathbf{Z}$ are orthogonal matrices. Meanwhile, $\mathbf{S}$ is the singular value matrix. A soft-thresholding operator is applied to singular values:

$$ \tilde{\mathbf{S}} = \text{diag}(\max(\text{diag}(\mathbf{S}) - \frac{\lambda_1}{\mu}, 0)), \quad (15) $$

where $\text{diag}(\cdot)$ returns a column vector consisting of the main diagonal elements if the input variable is a matrix and returns a diagonal matrix with the elements on the main diagonal if the input variable is a vector. $\mathbf{V}_i^{k+1}$ can then be expressed as

$$ \mathbf{V}_i^{k+1} = \tilde{\mathbf{U}} \tilde{\mathbf{S}} \mathbf{Z}^T. \quad (16) $$

Update to $\mathbf{E}$. The optimization problem for $\mathbf{E}_j$ is

$$ \begin{aligned} \mathbf{E}_j^{k+1} &= \arg \min_{\mathbf{E}_j^k} \lambda_2 \|\mathbf{E}_j^k\|_* + \frac{\mu}{2} \|\mathcal{B}_{:,;,j}^{k+1} - \mathbf{E}_j^k - \mathbf{H}_j^k\|_F^2 \\ &= \frac{1}{2} \|\tilde{\mathbf{E}}_j - \mathbf{E}_j^k\|_F^2 + \frac{\lambda_2}{\mu} \|\mathbf{E}_j^k\|_*; \end{aligned} \quad (17) $$

where $\tilde{\mathbf{E}}_j = \mathcal{B}_{:,;,j}^{k+1} - \mathbf{H}_j^k \in \mathbb{R}^{n_{row} \times n_{col}}$. Sub-problem (17) can be solved via SVD of $\tilde{\mathbf{E}}_j$:

$$ [\mathbf{U}, \mathbf{S}, \mathbf{Z}] = \text{SVD}(\tilde{\mathbf{E}}_j), \quad (18) $$

where $\mathbf{S}$ is the singular value matrix. As for (15), a soft-thresholding operator is applied to the singular values:

$$ \tilde{\mathbf{S}} = \text{diag}(\max(\text{diag}(\mathbf{S}) - \frac{\lambda_2}{\mu}, 0)), \quad (19) $$

then $\mathbf{E}_j^{k+1}$ can be expressed as

$$ \mathbf{E}_j^{k+1} = \tilde{\mathbf{U}} \tilde{\mathbf{S}} \mathbf{Z}^T. \quad (20) $$

Update $\mathbf{D}_i$.

$$ \mathbf{D}_i^{k+1} = \mathbf{D}_i^k - (\mathcal{A}_{:,;,i}^{k+1} - \mathbf{V}_i^{k+1}). \quad (21) $$

Update $\mathbf{H}_j$.

$$ \mathbf{H}_j^{k+1} = \mathbf{H}_j^k - (\mathcal{B}_{:,;,j}^{k+1} - \mathbf{E}_j^{k+1}). \quad (22) $$

Update $\mathbf{G}$.

$$ \mathbf{G}^{k+1} = \mathbf{G}^k - (\sum_{i=1}^{R} \mathcal{A}_{:,;,i}^{k+1} - \mathbf{P}). \quad (23) $$

### Algorithm 1: The Proposed LR-NTF Algorithm

**Input:** $\mathcal{Y}$; Matrix of Endmembers: $\mathbf{E}$

**Output:** A cube of abundance maps: $\mathcal{A}$

1. for $k=1; k < \text{Iter}; k++ \text{ do}$

2. $\mathcal{A}_{:,;,i}^{k+1} = (\sum_{b=1}^{L} c_{b_i} c_{b_i}^T + 2\mu \mathbf{I})^{-1} (\sum_{b=1}^{L} \mathcal{O}_{:,;,b} c_{b_i}^T + \mu (\mathbf{V}_i^k + \mathbf{D}_i^k + \mathbf{P} + \mathbf{G}^k - \tilde{\mathbf{A}}));$

3. $\mathcal{A} = \text{abs}(\mathcal{A})$;

4. $\mathcal{B}_{:,;,j}^{k+1} = (\sum_{b=1}^{L} m_{b_j} m_{b_j}^T + \mu \mathbf{I})^{-1} (\sum_{b=1}^{L} \mathcal{K}_{:,;,b} m_{b_j}^T + \mu (\mathbf{E}_j^k + \mathbf{H}_j^k)).$ If any element of $\mathcal{B}$ exceeds that of $\mathcal{A}^*$, it is replaced with that of $\mathcal{A}^*$;

5. $\mathcal{B} = \text{abs}(\mathcal{B})$;

6. $\mathbf{V}_i^{k+1} = \tilde{\mathbf{U}} \tilde{\mathbf{S}} \mathbf{Z}^T$;

7. $\mathbf{E}_j^{k+1} = \tilde{\mathbf{U}} \tilde{\mathbf{S}} \mathbf{Z}^T$;

8. $\mathbf{D}_i^{k+1} = \mathbf{D}_i^k - (\mathcal{A}_{:,;,i}^{k+1} - \mathbf{V}_i^{k+1});$

9. $\mathbf{H}_j^{k+1} = \mathbf{H}_j^k - (\mathcal{B}_{:,;,j}^{k+1} - \mathbf{E}_j^{k+1});$

10. $\mathbf{G}^{k+1} = \mathbf{G}^k - (\sum_{i=1}^{R} \mathcal{A}_{:,;,i}^{k+1} - \mathbf{P});$

11. $k = k + 1;$

12. end

13. return result
### III. EXPERIMENTS ON SYNTHETIC DATA

In this section, we compare the performance of the proposed algorithm LR-NTF and other algorithms including the GDA and semi-NMF. The GDA is considered to be the benchmark for solving the generalized bilinear model, and semi-NMF solves the unmixing problem in matrix form.

Two widely used metrics, namely the root-mean-square error (RMSE) of abundances and the image reconstruction error (RE) were adopted for evaluating the unmixing methods. The RMSE quantifies the difference between the estimated abundances $\mathcal{A}$ and the true abundances $\mathcal{A}$ as follows:

$$RMSE = \sqrt{\frac{1}{R \times N} \|\mathcal{A} - \hat{\mathcal{A}}\|_F^2} \quad (24)$$

The RE measures the difference between the observations $\mathcal{Y}$ and their reconstructions $\hat{\mathcal{Y}}$ as follows:

$$RE = \sqrt{\frac{1}{N \times L} \|\mathcal{Y} - \hat{\mathcal{Y}}\|_F^2} \quad (25)$$

#### A. Data Generation

In this study, the simulated data were synthesized in a similar way to the data in [33], and [30]. The synthesis was carried out as follows.

1) Six spectral reflection signals with 224 spectral bands ranging from 0.38 to 2.5 µm were chosen from the United States Geological Survey (USGS) digital spectral library<fnref n="1" />. Specifically, these were Carnallite, Ammonio-jarosite, Almandine, Brucite, Axinite, and Chlonte.

2) We generated an image of size $s^2 \times s^2 \times L$: this could be divided into small blocks of size $s \times s \times L$.

3) A randomly selected endmember was assigned to each block, and a $k \times k$ low-pass filter was then applied to generate abundance maps of size $s^2 \times s^2 \times R$ that contained mixed pixels while still satisfying the ANC and ASC constraints.

4) To simulate a non-pure pixel scenario, we detected the abundances whose values were higher than 0.8 and replaced them with the average fraction for all the endmembers.

5) The abundance and endmember information for the image were obtained in steps 1-4. Clean HSIs were generated based on two kinds of nonlinear mixing models, namely the generalized bilinear model and the polynomial post nonlinear model. The interaction coefficients in the GBM were set randomly, and the interaction coefficients in the PPNM were set to 0.25.

6) To evaluate the robustness to additive noise, zero-mean Gaussian white noise was added to the clean data. A series of noisy images with signal-to-noise ratios (SNRs) = {15, 20, 30, 40} dB was generated.

#### B. Algorithm Evaluation

After obtaining the synthetic data, we compared different unmixing methods, including the proposed method, the GDA,

and semi-NMF. Meanwhile, the proposed method was compared with the minimum volume nonnegative tensor factorization (MV-NTF) [30] separately on the three nonlinear datasets to see the impact of low-rank representation of interaction abundance maps. We first tested the effect of different SNRs before comparing the methods. The parameter settings for the tested methods are described in the following.

1) **Parameter Setting:** In this experiment, we analysed the parameters used in all the algorithms. For the GDA, which was the benchmark method, the tolerance for stopping the iterations was set to $1 \times 10^{-6}$. The parameters in the semi-NMF and MV-NTF were hand-tuned to get the optimal performance. The proposed LR-NTF has an abundance low-rank regularization parameter $\lambda_1$, a nonlinear iteration abundance low-rank regularization parameter $\lambda_2$, and a balance parameter $\mu$ that controls the convergence rate of the algorithm. We searched for the optimal value of $\lambda_1$ and $\lambda_2$ using a grid method with $\lambda_1 \in \{0.1, 0.2, ..., 1\}$ and $\lambda_2 \in \{0.01, 0.02, ..., 0.1\}$. To test the impact of parameters $\lambda_1$ and $\lambda_2$ on the unmixing, we generated a synthetic dataset with parameters $k = 9$, $s = 6$, and SNR = 30 dB. In our experiments, similar results were obtained for certain ranges of parameters: i.e. $\lambda_1 \in \{0.1\}$ and $\lambda_2 \in \{0.06, 0.07, 0.08, 0.09\}$. Hence, we set the parameters $\lambda_1$ and $\lambda_2$ to within this range in the subsequent experiments.

2) **Comparison between Methods Using Different Gaussian Noise Levels:** In this experiment, three images (of size $100 \times 100 \times 224$) were generated using two types of model (GBM and PPNM): image1 and image2, which are described in Table I, were generated using the GBM and PPNM, respectively. image3 was a mixture of image1 and image2, with half the pixels being generated by the GBM and the rest generated by the PPNM [16]. The first image, generated based on GBM, is used to compare the proposed algorithm to the GDA and semi-NMF method, whereas image2 and image3 were used to test the robustness of the proposed algorithm. Meanwhile, white Gaussian noise was added in the simulated images: the signal-to-noise ratio (SNR) was set to 15 dB, 20 dB, 30 dB, and 40 dB.

For all of the methods, the abundance matrix estimated by the FCLS algorithm [34] was used for the abundance initialization. Furthermore, as this was a supervised unmixing task, the endmember matrix was supposed to be known beforehand. Since the quality of the endmember extraction affects the results of the abundance inversion, for all the scenes shown in Table I, the true endmember information was used. The interaction endmember signals were obtained from the corresponding Hadamard product of two true endmembers. The parameters of LR-NTF were set as follows: $\lambda_1 = 0.1$, $\lambda_2 = 0.07$, and $\mu = 8 \times 10^{-3}$. Meanwhile, the maximum number of program iterations was set to 1000. We remark that theoretically, under the same constraints, the optimal solution of LR-NTF (whose objective function is a combination of the reconstruction error and two regularization terms) is expected to yield a higher reconstruction error than the solution of the Semi-NMF (whose objective function is the single reconstruction error). However, the exact optimal solutions for both above two problems are hard to obtain in practice. Semi-NMF and LR-NTF obtain approximated optimal solutions. And the

<footnote n="1">¹https://www.usgs.gov/labs/spec-lab</footnote>
<table><caption>TABLE I
THE EVALUATE RESULTS OF THE PROPOSED METHOD AND STATE-OF-THE-ART METHODS</caption><thead><tr><th>Scenario</th><th>SNR</th><th>Metric</th><th>FCLS</th><th>GDA</th><th>Semi-NMF</th><th>LR-NTF (ours)</th></tr></thead><tbody><tr><td rowspan="12">image1</td><td rowspan="3">15</td><td>RMSE</td><td>0.0746</td><td>0.0627</td><td>0.0608</td><td><strong>0.0437</strong></td></tr><tr><td>RE</td><td>0.0983</td><td>0.0975</td><td>0.0959</td><td><strong>0.0954</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1816</td><td>18</td><td>1022</td></tr><tr><td rowspan="3">20</td><td>RMSE</td><td>0.0680</td><td>0.0534</td><td>0.0473</td><td><strong>0.0253</strong></td></tr><tr><td>RE</td><td>0.0582</td><td>0.0570</td><td>0.0542</td><td><strong>0.0537</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1923</td><td>25</td><td>1246</td></tr><tr><td rowspan="3">30</td><td>RMSE</td><td>0.0646</td><td>0.0485</td><td>0.0322</td><td><strong>0.0146</strong></td></tr><tr><td>RE</td><td>0.0278</td><td>0.0252</td><td>0.0178</td><td><strong>0.0170</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1930</td><td>48</td><td>1330</td></tr><tr><td rowspan="3">40</td><td>RMSE</td><td>0.0641</td><td>0.0479</td><td>0.0267</td><td><strong>0.0141</strong></td></tr><tr><td>RE</td><td>0.0227</td><td>0.0194</td><td>0.0071</td><td><strong>0.0056</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1911</td><td>66</td><td>706</td></tr><tr><td rowspan="12">image2</td><td rowspan="3">15</td><td>RMSE</td><td>0.1050</td><td>0.0881</td><td>0.0755</td><td><strong>0.0453</strong></td></tr><tr><td>RE</td><td>0.1081</td><td>0.1066</td><td>0.1015</td><td><strong>0.1007</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1630</td><td>29</td><td>833</td></tr><tr><td rowspan="3">20</td><td>RMSE</td><td>0.1011</td><td>0.0815</td><td>0.0627</td><td><strong>0.0305</strong></td></tr><tr><td>RE</td><td>0.0686</td><td>0.0663</td><td>0.0579</td><td><strong>0.0568</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1779</td><td>40</td><td>906</td></tr><tr><td rowspan="3">30</td><td>RMSE</td><td>0.0993</td><td>0.0788</td><td>0.0511</td><td><strong>0.0233</strong></td></tr><tr><td>RE</td><td>0.0426</td><td>0.0386</td><td>0.0210</td><td><strong>0.0185</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1841</td><td>62</td><td>585</td></tr><tr><td rowspan="3">40</td><td>RMSE</td><td>0.0991</td><td>0.0785</td><td>0.0489</td><td><strong>0.0224</strong></td></tr><tr><td>RE</td><td>0.0390</td><td>0.0347</td><td>0.0121</td><td><strong>0.0075</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1817</td><td>66</td><td>1349</td></tr><tr><td rowspan="12">image3</td><td rowspan="3">15</td><td>RMSE</td><td>0.0910</td><td>0.0769</td><td>0.0689</td><td><strong>0.0444</strong></td></tr><tr><td>RE</td><td>0.1031</td><td>0.1020</td><td>0.0989</td><td><strong>0.0980</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1735</td><td>24</td><td>880</td></tr><tr><td rowspan="3">20</td><td>RMSE</td><td>0.0854</td><td>0.0685</td><td>0.0550</td><td><strong>0.0286</strong></td></tr><tr><td>RE</td><td>0.0634</td><td>0.0616</td><td>0.0560</td><td><strong>0.0552</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1895</td><td>34</td><td>1082</td></tr><tr><td rowspan="3">30</td><td>RMSE</td><td>0.0832</td><td>0.0653</td><td>0.0417</td><td><strong>0.0198</strong></td></tr><tr><td>RE</td><td>0.0356</td><td>0.0324</td><td>0.0194</td><td><strong>0.0178</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1878</td><td>58</td><td>677</td></tr><tr><td rowspan="3">40</td><td>RMSE</td><td>0.0830</td><td>0.0651</td><td>0.0387</td><td><strong>0.0186</strong></td></tr><tr><td>RE</td><td>0.0315</td><td>0.0279</td><td>0.0099</td><td><strong>0.0067</strong></td></tr><tr><td>Time(s)</td><td><strong>3</strong></td><td>1893</td><td>66</td><td>667</td></tr></tbody></table>

<table><caption>TABLE II
THE EVALUATE RESULTS OF THE PROPOSED METHOD AND THE MV-NTF</caption><thead><tr><th>Scenario</th><th>SNR</th><th>Metric</th><th>MV-NTF</th><th>LR-NTF-MVNTF (ours)</th></tr></thead><tbody><tr><td rowspan="12">image1</td><td rowspan="3">15</td><td>RMSE</td><td>0.0831</td><td><strong>0.0816</strong></td></tr><tr><td>RE</td><td>0.0984</td><td><strong>0.0951</strong></td></tr><tr><td>Time(s)</td><td>1299</td><td><strong>489</strong></td></tr><tr><td rowspan="3">20</td><td>RMSE</td><td>0.0693</td><td><strong>0.0669</strong></td></tr><tr><td>RE</td><td>0.0575</td><td><strong>0.0534</strong></td></tr><tr><td>Time(s)</td><td>1301</td><td><strong>506</strong></td></tr><tr><td rowspan="3">30</td><td>RMSE</td><td>0.0713</td><td><strong>0.0643</strong></td></tr><tr><td>RE</td><td>0.0262</td><td><strong>0.0170</strong></td></tr><tr><td>Time(s)</td><td>1300</td><td><strong>708</strong></td></tr><tr><td rowspan="3">40</td><td>RMSE</td><td>0.0719</td><td><strong>0.0649</strong></td></tr><tr><td>RE</td><td>0.0198</td><td><strong>0.0057</strong></td></tr><tr><td>Time(s)</td><td>1293</td><td><strong>875</strong></td></tr><tr><td rowspan="12">image2</td><td rowspan="3">15</td><td>RMSE</td><td><strong>0.0752</strong></td><td>0.0791</td></tr><tr><td>RE</td><td>0.1029</td><td><strong>0.1001</strong></td></tr><tr><td>Time(s)</td><td>1297</td><td><strong>494</strong></td></tr><tr><td rowspan="3">20</td><td>RMSE</td><td>0.0713</td><td><strong>0.0694</strong></td></tr><tr><td>RE</td><td>0.0587</td><td><strong>0.0563</strong></td></tr><tr><td>Time(s)</td><td>1298</td><td><strong>609</strong></td></tr><tr><td rowspan="3">30</td><td>RMSE</td><td>0.0662</td><td><strong>0.0589</strong></td></tr><tr><td>RE</td><td>0.0221</td><td><strong>0.0179</strong></td></tr><tr><td>Time(s)</td><td>1304</td><td><strong>680</strong></td></tr><tr><td rowspan="3">40</td><td>RMSE</td><td>0.0766</td><td><strong>0.0707</strong></td></tr><tr><td>RE</td><td>0.0155</td><td><strong>0.0060</strong></td></tr><tr><td>Time(s)</td><td>1304</td><td><strong>652</strong></td></tr><tr><td rowspan="12">image3</td><td rowspan="3">15</td><td>RMSE</td><td><strong>0.0800</strong></td><td>0.0847</td></tr><tr><td>RE</td><td>0.1013</td><td><strong>0.0975</strong></td></tr><tr><td>Time(s)</td><td>1299</td><td><strong>465</strong></td></tr><tr><td rowspan="3">20</td><td>RMSE</td><td>0.0683</td><td><strong>0.0657</strong></td></tr><tr><td>RE</td><td>0.0592</td><td><strong>0.0548</strong></td></tr><tr><td>Time(s)</td><td>1297</td><td><strong>536</strong></td></tr><tr><td rowspan="3">30</td><td>RMSE</td><td>0.0722</td><td><strong>0.0650</strong></td></tr><tr><td>RE</td><td>0.0277</td><td><strong>0.0174</strong></td></tr><tr><td>Time(s)</td><td>1297</td><td><strong>755</strong></td></tr><tr><td rowspan="3">40</td><td>RMSE</td><td>0.0676</td><td><strong>0.0597</strong></td></tr><tr><td>RE</td><td>0.0217</td><td><strong>0.0058</strong></td></tr><tr><td>Time(s)</td><td>1305</td><td><strong>957</strong></td></tr></tbody></table>

applied to image1, but also for image2 and image3.

feasibility of these solutions is not guaranteed. So that it is possible for the solution of LR-NTF to be better in the aspect of the reconstruction error.

It can be seen from Fig. 3 that, for the different levels of noise in image1, the proposed method uniformly yields the best performance in terms of RMSE and RE. Of the three nonlinear unmixing algorithms, the GDA produces the worst results since it carries out the unmixing pixel by pixel without considering spatial correlation or sparsity of abundances. The semi-NMF and the proposed LR-NTF use nonnegative matrix factorization and tensor factorization, respectively, to process a whole image. These two kinds of factorization exploit the spatial correlation of HSIs. Furthermore, the proposed method takes advantage of the low rank of abundance maps, which boosts its unmixing performance.

To evaluate the robustness of the proposed algorithm to the different types of spectral mixing, we generated three images: image1, which was based on the GBM, image2 based on the PPNM, and image3, for which half the pixels were based on the GBM and the other half were based on the PPNM. It can be seen from Table I that, although it was based on the GBM model, the proposed LR-NTF performed best not only when

<figure><img src="image_6.png" /><figcaption>Fig. 3. Unmixing performance in terms of RMSE (a) and RE (b) in the simulated image1 with different SNRs</figcaption></figure>

To evaluate the impact of low-rank representations of interaction abundance maps in nonlinear mixing images, we conducted a comparison between our method and the MV-NTF method [30], which unmixes image based on LMM model by using low-rank representations of abundance maps. MV-NTF method estimates abundances and endmembers jointly,
<figure><img src="image_7.png" /><figcaption>Fig. 4. Iterations of RMSE (a) and RE (b) in LR-NTF with the simulated image1</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Iterations of RMSE (a) and RE (b) in LR-NTF with the simulated image2</figcaption></figure>

whereas the proposed LR-NTF only estimates abundances. For a fair comparison, a modified version of LR-NTF method, called 'LR-NTF-MVNTF', was implemented using endmembers estimated by MV-NTF, instead of being estimated by VCA. Results in Table II show that our method, LR-NTF-MVNTF, can obtain the best results in most cases with gains increasing as the SNR increases.

Fig. 4, Fig. 5,and Fig. 6 show the iterations of RMSE and RE with the image1 (30 dB), image2 (30 dB), image3 (30 dB),

<figure><img src="image_9.png" /><figcaption>Fig. 6. Iterations of RMSE (a) and RE (b) in LR-NTF with the simulated image3</figcaption></figure>

respectively [16]. It is shown that the proposed algorithm is an excellent algorithm with good convergence.

As can be seen in Table I, the FCLS obtains the least run time, and the GDA obtains the longest running time. Although the proposed method is based on the NTF, it's also faster than GDA and MV-NTF.

## IV. EXPERIMENTS WITH REAL DATASET

In this section, We describe tests where the proposed tensor-based nonlinear unmixing algorithm was applied to four real images. Two of the images were acquired by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS), and the others were collected by the Hyperspectral Digital Imagery Collection Experiment (HYDICE) sensor. Due to the lack of the ground truth for the abundances, the RE in (25) and average of spectral angle mapper (aSAM) were used to measure the performance of both the proposed method and the state-of-the-art methods. The aSAM metric can qualify the average spectral angle mapping of the reconstructed jth spectral vector $\hat{y}_j$ and observed jth spectral vector $y_j$:

$$ aSAM = \frac{1}{N} \sum_{j=1}^{N} \arccos \left( \frac{\mathbf{y}_j^T \cdot \hat{\mathbf{y}}_j}{\|\mathbf{y}_j\| \|\hat{\mathbf{y}}_j\|} \right). \quad (26) $$

Note that RE and aSAM are reference metrics, not direct metrics measuring the quality of estimated abundances.

### A. Cuprite

The first of the real datasets was the widely used hyperspectral image for testing unmixing methods [35]–[38], which was acquired by the AVIRIS sensor over the Cuprite mining region, NV, USA. The image has 224 spectral bands in the range 0.4 ~ 2.5µm. After removing the bands with a low SNR and the water absorption bands, 188 bands remained. A sub-image with 250 × 191 pixels was chosen for use in our experiments: refer to Fig. 7(a) for details. The spatial distribution of minerals can be inferred from a reference mineral map (shown in Fig. 7(b)) created by a classifier based on domain-specific knowledge.

The sub-image that was used has been extensively studied in [35]–[38]. We set the number of endmembers, R, to 10, the same as in [35], and [39]. Meanwhile, VCA [31] was used to extract the endmember signatures, which were Andradite, Sphene, Muscovite, Montmorillonite CM20, Kaolinite, Alunite, Nontronite, Montmorillonite Na, Buddingtonite, and Chalcedony.

The ASC regularization term in the semi-NMF and MV-NTF were set to 0.1 and 0.6, respectively. For the proposed LR-NTF method, the low-rank abundance parameter and the nonlinear interaction regularization terms were set to 0.1 and 0.07, respectively. Furthermore, the penalty parameter $\mu = 1 \times 10^{-4}$, and the abundance maps were initialized using the FCLS algorithm. The tolerance level, at which the iterations were stopped, was set to $1 \times 10^{-6}$ for GDA, semi-NMF, and LR-NTF. The number of iterations of the proposed LR-NTF was 260 for the Cuprite image (see Fig. 8(b)).

Table III shows the RE, aSAM, and time cost for all the algorithms. It can be seen that, as it considers the spatial and
<figure><img src="image_10.png" /><figcaption>Fig. 7. HSIs used in the experiments: (a) sub-image of Cuprite provided by USGS in 1995, (b) reference map of Cuprite, (c) sub-image of San Diego Airport data, (d) sub-image1 of Washington Dc Mall, (e) sub-image2 of Washington Dc Mall</figcaption></figure>

<table><caption>TABLE III
EVALUATION RESULTS IN CUPRITE WITH RE, ASAM AND TIME COST (S).</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Metric</th><th colspan="4">Using endmembers extracted by VCA</th><th colspan="2">Using endmembers extracted by MVNTF</th></tr><tr><th>FCLS</th><th>GDA</th><th>Semi-NMF</th><th>LR-NTF (ours)</th><th>MV-NTF</th><th>LR-NTF-MVNTF (ours)</th></tr></thead><tbody><tr><td rowspan="3">Cuprite</td><td>RE</td><td>0.0080</td><td>0.0080</td><td>0.0061</td><td><strong>0.0053</strong></td><td>0.0104</td><td><strong>0.0038</strong></td></tr><tr><td>aSAM</td><td>0.0191</td><td>0.0191</td><td>0.0165</td><td><strong>0.0145</strong></td><td>0.0134</td><td><strong>0.0107</strong></td></tr><tr><td>Time</td><td><strong>35</strong></td><td>176</td><td>93</td><td>397</td><td>4643</td><td>187</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Fig. 8. Image reconstruction errors (REs) and the convergence of LR-NTF in Cuprite data</figcaption></figure>

low-rank information in the HSIs simultaneously, the error produced by the proposed LR-NTF is the lowest. Meanwhile, the proposed LR-NTF using endmembers extracted by MV-NTF obtains lower RE and aSAM than the LR-NTF using endmembers extracted by VCA. Fig. 9 and Fig. 10 show the estimated abundance maps and distributions of the REs, respectively. The significant differences of our results compared to other methods in Fig. 9 are marked in carmine circles. The bright areas in Fig. 10 indicate large errors in the reconstructed images. The errors for the FCLS are the worst because this method only considers linear mixtures of the minerals. The semi-NMF performs better than GDA because the GDA is a pixel-based algorithm that does not take any spatial information into consideration. As can be seen from Fig. 10, the proposed algorithm performs significantly better than the semi-NMF as it takes into account spatial information from the region of interest by using the tensor-based form and low-rank distribution of materials.

## B. San Diego Airport

The second real dataset used in the experiment was an airborne hyperspectral image of San Diego Airport that was acquired by the AVIRIS sensor. This image has been widely used in hyperspectral target detection [40], [41]. The image contains 400 × 400 pixels and 224 spectral bands. As in [42], a sub-image with 160 × 140 pixels was chosen for our experiments, and is shown in Fig. 7(c). In order to avoid the errors caused by water vapor absorption and low SNRs, bands 1-6, 33-35, 97, 107-113, 153-166, and 221-224 were removed. The sub-image included the following five main materials: Roof, Grass, Ground and Road, Tree, and Other (see [42] for more details). We, therefore, set $R = 5$ in the experiments; the VCA was used to extract the endmembers.

The optimal parameter values were used in the comparison experiments. The ANC balance parameters in the semi-NMF and MV-NTF were set to 0.1 and 0.6, respectively. For the proposed algorithm, the low-rank abundance regularization parameter was set to 0.1, the low-rank nonlinear interaction regularization parameter was set to 0.07, and the penalty parameter to $1 \times 10^{-4}$. Furthermore, the FCLS was used to

initialize the abundance map, and the maximum number of iterations was set to 1000. The tolerance used for stopping the iterations of a solver was set to $1 \times 10^{-6}$ for GDA, semi-NMF, and LR-NTF. Fig. 11(b) shows a plot of RE against the number of iterations for the San Diego Airport data. Table IV shows the RE and aSAM for all the algorithms. It can be seen that, as it considers the spatial and low-rank information in the HSIs simultaneously, the error produced by the proposed LR-NTF is the lowest. A comparison between results of LR-NTF (using endmembers estimated by MV-NTF) and LR-NTF-MVNTF using endmembers estimated by VCA) indicates that the performance of the proposed LR-NTF can be improved by using a higher quality of endmembers.

It can be seen from Fig. 12 that the proposed algorithm successfully distinguishes between the Roof and the Ground&Road classes marked by carmine circles. Fig. 11(a) shows the RE values for all of the algorithms, and it can be seen that the proposed algorithm produces the best results. Further details are given in Fig. 13, where the RE distribution is shown visually. The very bright areas in the image correspond to large REs: it is clear from this image that the error map obtained using the proposed algorithm shows smaller error than the ones produced by the other algorithms.

## C. Washington DC Mall

The third real dataset used in the experiment was the Washington DC Mall scene, which was acquired by the HYDICE sensor over Washington DC, USA. The full image contains 1208 × 307 pixels, with 210 spectral bands ranging from 0.4 to 2.5 µm. The spatial resolution is 3 m. After removing the atmospheric absorption and low-SNR bands (bands 103–106, 138–148, and 207–210), 191 bands remained. This image contains seven main materials: Roof, Grass, Road, Trail, Water, Shadow, and Tree. Two sub-images were chosen from the original image to test the proposed algorithms.

The first sub-image consisting of 256x256 pixels, called sub-DC1, was clipped from the Washington DC Mall data (see Fig. 7(d)). Hysime [43] and the VCA were used to estimate the number of endmembers and the endmember matrix, respectively. The extracted endmembers were called Roof1, Roof2, Grass, Road, Tree, and Trail. The optimal parameter values for all of the algorithms were set as follows. The ASC parameter was set to 0.1 in the semi-NMF; for the proposed LR-NTF method, the low-rank regularization parameter was set to 0.1 in the abundance map and 0.07 in the nonlinear interaction map. The LR-NTF penalty parameter was set to $1 \times 10^{-4}$, and FCLS was used to initialize the abundance maps. The tolerance level used to stop the iterations of a solver was set to $1 \times 10^{-6}$ in GDA, semi-NMF, and LR-NTF. Fig. 14 shows the RE for all the algorithms together with the RE plotted against the number of iterations of the proposed algorithm. Fig. 15 depicts the estimated abundance maps for the proposed algorithm as well as the other algorithms. A detailed analysis shows that the FCLS produces poor estimates of the abundance maps since the LMM fails to model the complex information in the scene. Meanwhile, the semi-NMF performs better than the GDA as it takes into account some
<figure><img src="image_12.png" /><figcaption>Fig. 9. Estimated abundance maps of the Cuprite image</figcaption></figure>

<table><caption>TABLE IV
EVALUATION RESULTS IN SAN DIEGO AIRPORT WITH RE, ASAM AND TIME COST (S).</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Metric</th><th colspan="4">Using endmembers extracted by VCA</th><th colspan="2">Using endmembers extracted by MVNTF</th></tr><tr><th>FCLS</th><th>GDA</th><th>Semi-NMF</th><th>LR-NTF (ours)</th><th>MV-NTF</th><th>LR-NTF-MVNTF (ours)</th></tr></thead><tbody><tr><td rowspan="3">San Diego Airport</td><td>RE</td><td>0.0165</td><td>0.0164</td><td>0.0150</td><td><strong>0.0132</strong></td><td>0.0184</td><td><strong>0.0110</strong></td></tr><tr><td>aSAM</td><td>0.0596</td><td>0.0594</td><td>0.0542</td><td><strong>0.0455</strong></td><td>0.0588</td><td><strong>0.0412</strong></td></tr><tr><td>Time</td><td><strong>4</strong></td><td>193</td><td>50</td><td>62</td><td>4946</td><td>188</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Fig. 10. RE distribution maps comparison between the proposed algorithm and state-of-the-art algorithms on the Cuprite</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. RE of four unmixing methods in (a) and iteration of RE in the LR-NTF in (b) with the San Diego Airport data</figcaption></figure>

of the spatial information by using the hyperspectral matrix. Furthermore, by considering the low-rank representation and the third-order tensor, the proposed method produces effective results. In order to study the differences in the error produced by all the algorithms, the RE distributions are shown visually
<table><caption>TABLE V
EVALUATION RESULTS IN WASHINGTON DC MALL WITH RE, ASAM AND TIME COST (S).</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Metric</th><th colspan="4">Using endmembers extracted by VCA</th><th colspan="2">Using endmembers extracted by MVNTF</th></tr><tr><th>FCLS</th><th>GDA</th><th>Semi-NMF</th><th>LR-NTF (ours)</th><th>MV-NTF</th><th>LR-NTF-MVNTF (ours)</th></tr></thead><tbody><tr><td rowspan="3">SubDc1</td><td>RE</td><td>0.0156</td><td>0.0154</td><td>0.0120</td><td><strong>0.0099</strong></td><td>0.0194</td><td><strong>0.0096</strong></td></tr><tr><td>aSAM</td><td>0.1020</td><td>0.1015</td><td>0.0837</td><td><strong>0.0623</strong></td><td>0.0980</td><td><strong>0.0636</strong></td></tr><tr><td>Time</td><td><strong>15</strong></td><td>765</td><td>45</td><td>381</td><td>11664</td><td>1370</td></tr><tr><td rowspan="3">SubDc2</td><td>RE</td><td>0.0094</td><td>0.0093</td><td>0.0080</td><td><strong>0.0078</strong></td><td>0.0071</td><td><strong>0.0027</strong></td></tr><tr><td>aSAM</td><td>0.0264</td><td>0.0264</td><td>0.0246</td><td><strong>0.0239</strong></td><td>0.0152</td><td><strong>0.0076</strong></td></tr><tr><td>Time</td><td><strong>7</strong></td><td>636</td><td>20</td><td>218</td><td>3355</td><td>771</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Estimated abundance maps comparison between the proposed algorithm and state-of-the-art algorithms on the San Diego Airport</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. RE distribution maps comparison between the proposed algorithm and state-of-the-art algorithms on the San Diego Airport</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. RE of four unmixing methods in (a) and iteration of RE in the LR-NTF in (b) with the sub-image1 of Washington DC Mall data</figcaption></figure>

in Fig. 16. The results in Table V illustrates that the proposed algorithm achieves a smaller error than other methods. After using the endmembers extracted by MV-NTF, the proposed LR-NTF obtains the more reasonable abundance estimation results than that by VCA.

The second sub-image, called sub-DC2, consisted of

180×160 pixels and was also clipped from the Washington DC Mall data (see Fig. 7(e)). This image has been studied in [42], and six materials were extracted: Tree, Trail, Roof, Water, Grass, and Road. Hence, the number of endmembers, R, was set as 6 for our experiment; the endmember matrix was extracted using the VCA, and the parameter values used for all the algorithms were the same as for the first sub-image. Fig. 17(a) shows that the proposed algorithm produces the lowest image reconstruction error and Fig. 18 shows the estimated abundance maps. By taking into account the low-rank characteristics of the abundance maps and the third-order representation of the HSI, the proposed algorithm produces the best estimation of the abundance maps. Fig. 19 shows the error distribution for the whole image. The FCLS performs worst since this method fails to dig the complex information from the scene. The semi-NMF performs better than the GDA as it uses a matrix-based technique instead of a pixel-based unmixing method. Furthermore, the proposed algorithm produces a smaller error map than the FCLS, GDA, or semi-NMF.

In the experiments with real images, we have tested the proposed method on the images with diverse spatial resolutions. Our method achieves good results in high spatial resolution images, including Washington DC Mall with spatial resolution 3m and San Diego Airport image with spatial resolution 3.5m. The spatial resolution of Cuprite image (20m) is relatively low, but our method still achieves good results. A critical prior knowledge the proposed method used is the low-rankness of abundance maps and interaction abundance maps. The low-rankness exists because of the high spatial correlation
<figure><img src="image_18.png" /><figcaption>Fig. 15. Estimated abundance maps comparison between the proposed algorithm and state-of-the-art algorithms on sub-image1 of Washington DC Mall data</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. RE distribution maps comparison between the proposed algorithm and state-of-the-art algorithms on sub-image1 of Washington DC Mall data</figcaption></figure>

of materials/endmembers, and does not depend on the spatial resolution of the image. Consequently, our proposed method could be used for images with diverse spatial resolutions, even

lower spatial resolution than the ones discussed in the real data experiments.

## V. CONCLUSION

In this paper, we proposed a nonnegative tensor factorization-based nonlinear hyperspectral unmixing method. By taking full advantage of the low rank of abundance maps, we imposed the nuclear norm on abundance maps and non-linear interaction maps. In order to evaluate the effectiveness of the proposed algorithm, several kinds of synthetic datasets and four real hyperspectral datasets were tested. Our method exploits the low-rank of abundance maps, and it was shown
<figure><img src="image_20.png" /><figcaption>Fig. 17. RE of four unmixing methods in (a) and iteration of RE in the LR-NTF in (b) with the sub-image2 of Washington DC Mall data</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. Estimated abundance maps comparison between the proposed algorithm and state-of-the-art algorithms on sub-image2 of Washington DC Mall data</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. RE distribution maps comparison between the proposed algorithm and state-of-the-art algorithms on sub-image2 of Washington DC Mall data</figcaption></figure>

<table><thead><tr><td>Abk.</td><td>Bedeutung</td></tr></thead><tbody><tr><td>aSAM</td><td>average of spectral angle mapper</td></tr><tr><td>BMMs</td><td>bilinear mixing models</td></tr><tr><td>GBM</td><td>generalized bilinear model</td></tr><tr><td>GDA</td><td>gradient descent algorithm</td></tr><tr><td>HNU</td><td>hyperspectral nonlinear unmixing</td></tr><tr><td>HSCs</td><td>hyperspectral cameras</td></tr><tr><td>HSIs</td><td>hyperspectral images</td></tr><tr><td>HU</td><td>hyperspectral unmixing</td></tr><tr><td>LMM</td><td>linear mixing model</td></tr><tr><td>LQM</td><td>linearquadratic mixing model</td></tr><tr><td>MHPNMM</td><td>multiharmonic postnonlinear mixing model</td></tr><tr><td>MLM</td><td>multilinear mixing model</td></tr><tr><td>MTL</td><td>multi-task learning</td></tr><tr><td>NTF</td><td>nonnegative tensor factorization,</td></tr><tr><td>NLMMs</td><td>nonlinear mixing models</td></tr><tr><td>RE</td><td>reconstruction error</td></tr><tr><td>RMSE</td><td>root-mean-square error</td></tr><tr><td>semi-NMF</td><td>semi-nonnegative matrix factorization</td></tr><tr><td>VCA</td><td>vertex component analysis</td></tr></tbody></table>

that this can improve the nonlinear unmixing performance. Furthremore, we tested the proposed method in images with diverse number of endmembers (from 5 to 10), and obtained better results compared to other methods. However, with the growing of the number of endmembers and the size of the images, the difficulty of unmixing will also increase, such as the need of more computation and memory, which is our future research direction.

# APPENDIX A

## ACKNOWLEDGMENT

The authors would like to thank Professor Naoto Yokoya for providing the semi-NMF code for our comparison experiment. Professor Yuntao Qian provided the endmember data used in some of the experiments with synthetic data.

## REFERENCES

[1] N. Keshava and J. F. Mustard, “Spectral unmixing,” *IEEE Signal Processing Magazine*, vol. 19, no. 1, pp. 44–57, Jan. 2002.

[2] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches,” *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, vol. 5, no. 2, pp. 354–379, Apr. 2012.

[3] N. Dobigeon, J. Tourneret, C. Richard, J. C. M. Bermudez, S. McLaughlin, and A. O. Hero, “Nonlinear unmixing of hyperspectral images: Models and algorithms,” *IEEE Signal Processing Magazine*, vol. 31, no. 1, pp. 82–94, Jan. 2014.

[4] B. Zhang, L. Zhuang, L. Gao, W. Luo, Q. Ran, and Q. Du, “PSO-EM: A hyperspectral unmixing algorithm based on normal compositional model,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 52, no. 12, pp. 7782–7792, Dec. 2014.

[5] D. Hong, N. Yokoya, J. Chanussot, and X. X. Zhu, “An augmented linear mixing model to address spectral variability for hyperspectral unmixing,” *IEEE Transactions on Image Processing*, vol. 28, no. 4, pp. 1923–1938, Apr. 2019.

[6] J. Yao, D. Meng, Q. Zhao, W. Cao, and Z. Xu, “Nonconvex-sparsity and nonlocal-smoothness-based blind hyperspectral unmixing,” *IEEE Transactions on Image Processing*, vol. 28, no. 6, pp. 2991–3006, Jun. 2019.

[7] H. Yu, L. Gao, W. Liao, B. Zhang, L. Zhuang, M. Song, and J. Chanussot, “Global spatial and local spectral similarity-based manifold learning group sparse representation for hyperspectral imagery classification,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 58, no. 5, pp. 3043–3056, May. 2020.
[8] 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," *International Journal of Remote Sensing*, vol. 30, no. 11, pp. 2951-2962, 2009. [Online]. Available: https://doi.org/10.1080/01431160802558659

[9] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret, "Nonlinear unmixing of hyperspectral images using a generalized bilinear model," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 49, no. 11, pp. 4153-4162, Nov. 2011.

[10] I. Meganem, P. Déliot, X. Briottet, Y. Deville, and S. Hosseini, "Linear-quadratic mixing model for reflectances in urban environments," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 52, no. 1, pp. 544-558, Jan. 2014.

[11] R. Heylen and P. Scheunders, "A multilinear mixing model for nonlinear spectral unmixing," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 54, no. 1, pp. 240-251, Jan. 2016.

[12] A. Marinoni, J. Plaza, A. Plaza, and P. Gamba, "Estimating nonlinearities in p-linear hyperspectral mixtures," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 56, no. 11, pp. 6586-6595, Nov. 2018.

[13] M. Tang, B. Zhang, A. Marinoni, L. Gao, and P. Gamba, "Multiharmonic postnonlinear mixing model for hyperspectral nonlinear unmixing," *IEEE Geoscience and Remote Sensing Letters*, vol. 15, no. 11, pp. 1765-1769, Nov. 2018.

[14] Q. Qu, N. M. Nasrabadi, and T. D. Tran, "Abundance estimation for bilinear mixture models via joint sparse and low-rank representation," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 52, no. 7, pp. 4404-4423, Jul. 2014.

[15] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret, "Unmixing hyperspectral images using the generalized bilinear model," in *2011 IEEE International Geoscience and Remote Sensing Symposium*, Jul. 2011, pp. 1886-1889.

[16] N. Yokoya, J. Chanussot, and A. Iwasaki, "Nonlinear unmixing of hyperspectral data using semi-nonnegative matrix factorization," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 52, no. 2, pp. 1430-1437, Feb. 2014.

[17] C. Li, Y. Ma, J. Huang, X. Mei, C. Liu, and J. Ma, "Gbm-based unmixing of hyperspectral data using bound projected optimal gradient method," *IEEE Geoscience and Remote Sensing Letters*, vol. 13, no. 7, pp. 952-956, Jul. 2016.

[18] Y. Su, J. Li, H. Qi, P. Gamba, A. Plaza, and J. Plaza, "Multi-task learning with low-rank matrix factorization for hyperspectral nonlinear unmixing," in *IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium*, Jul. 2019, pp. 2127-2130.

[19] M. A. Veganzones, J. E. Cohen, R. Cabral Farias, J. Chanussot, and P. Comon, "Nonnegative tensor cp decomposition of hyperspectral data," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 54, no. 5, pp. 2577-2588, May. 2016.

[20] N. Renard, S. Bourennane, and J. Blanc-Talon, "Denoising and dimensionality reduction using multilinear tools for hyperspectral images," *IEEE Geoscience and Remote Sensing Letters*, vol. 5, no. 2, pp. 138-142, Apr. 2008.

[21] T. Imbiriba, R. A. Borsoi, and J. C. M. Bermudez, "A low-rank tensor regularization strategy for hyperspectral unmixing," in *2018 IEEE Statistical Signal Processing Workshop (SSP)*, Jun. 2018, pp. 373-377.

[22] F. Xiong, Y. Qian, J. Zhou, and Y. Y. Tang, "Hyperspectral unmixing via total variation regularized nonnegative tensor factorization," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 57, no. 4, pp. 2341-2357, Apr. 2019.

[23] P. V. Giampouras, K. E. Themelis, A. A. Rontogiannis, and K. D. Koutroumbas, "Simultaneously sparse and low-rank abundance matrix estimation for hyperspectral image unmixing," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 54, no. 8, pp. 4775-4789, Aug. 2016.

[24] Y. Xu, Z. Wu, J. Li, A. Plaza, and Z. Wei, "Anomaly detection in hyperspectral images based on low-rank and sparse representation," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 54, no. 4, pp. 1990-2000, Apr. 2016.

[25] J. Yang, Y. Zhao, J. C. Chan, and S. G. Kong, "Coupled sparse denoising and unmixing with low-rank constraint for hyperspectral image," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 54, no. 3, pp. 1818-1833, Mar. 2016.

[26] C. G. Tsinos, A. A. Rontogiannis, and K. Berberidis, "Distributed blind hyperspectral unmixing via joint sparsity and low-rank constrained non-negative matrix factorization," *IEEE Transactions on Computational Imaging*, vol. 3, no. 2, pp. 160-174, Jun. 2017.

[27] D. Hong and X. X. Zhu, "Sulora: Subspace unmixing with low-rank attribute embedding for hyperspectral data analysis," *IEEE Journal of Selected Topics in Signal Processing*, vol. 12, no. 6, pp. 1351-1363, Dec. 2018.

[28] F. Feng, B. Zhao, L. Tang, W. Wang, and S. Jia, "Robust low-rank abundance matrix estimation for hyperspectral unmixing," *The Journal of Engineering*, vol. 2019, no. 21, pp. 7406-7409, 2019.

[29] J. Huang, T. Huang, L. Deng, and X. Zhao, "Joint-sparse-blocks and low-rank representation for hyperspectral unmixing," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 57, no. 4, pp. 2419-2438, Apr. 2019.

[30] Y. Qian, F. Xiong, S. Zeng, J. Zhou, and Y. Y. Tang, "Matrix-vector nonnegative tensor factorization for blind unmixing of hyperspectral imagery," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 55, no. 3, pp. 1776-1792, Mar. 2017.

[31] J. M. P. Nascimento and J. M. B. Dias, "Vertex component analysis: a fast algorithm to unmix hyperspectral data," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 43, no. 4, pp. 898-910, Apr. 2005.

[32] J. Eckstein and D. P. Bertsekas, "On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators," *Mathematical Programming*, 1992.

[33] L. Miao and H. Qi, "Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 45, no. 3, pp. 765-777, Mar. 2007.

[34] D. C. Heinz and Chein-I-Chang, "Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 39, no. 3, pp. 529-545, Mar. 2001.

[35] Y. Qian, S. Jia, J. Zhou, and A. Robles-Kelly, "Hyperspectral unmixing via $l_1/2$ sparsity-constrained nonnegative matrix factorization," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 49, no. 11, pp. 4282-4297, Nov. 2011.

[36] Q. Wei, M. Chen, J. Tourneret, and S. Godsill, "Unsupervised nonlinear spectral unmixing based on a multilinear mixing model," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 55, no. 8, pp. 4534-4544, Aug. 2017.

[37] X. Wang, Y. Zhong, L. Zhang, and Y. Xu, "Spatial group sparsity regularized nonnegative matrix factorization for hyperspectral unmixing," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 55, no. 11, pp. 6287-6304, Nov. 2017.

[38] L. Zhuang, C. Lin, M. A. T. Figueiredo, and J. M. Bioucas-Dias, "Regularization parameter selection in minimum volume hyperspectral unmixing," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 57, no. 12, pp. 9858-9877, Dec. 2019.

[39] J. Li, J. M. Bioucas-Dias, A. Plaza, and L. Liu, "Robust collaborative nonnegative matrix factorization for hyperspectral unmixing," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 54, no. 10, pp. 6076-6090, Oct. 2016.

[40] A. Taghipour and H. Ghassemian, "Hyperspectral anomaly detection using attribute profiles," *IEEE Geoscience and Remote Sensing Letters*, vol. 14, no. 7, pp. 1136-1140, Jul. 2017.

[41] M. Vafadar and H. Ghassemian, "Anomaly detection of hyperspectral imagery using modified collaborative representation," *IEEE Geoscience and Remote Sensing Letters*, vol. 15, no. 4, pp. 577-581, Apr. 2018.

[42] F. Zhu, "Hyperspectral unmixing: Ground truth labeling, datasets, benchmark performances and survey," *arXiv preprint arXiv:1708.05125*, 2017.

[43] J. M. Bioucas-Dias and J. M. P. Nascimento, "Hyperspectral subspace identification," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 46, no. 8, pp. 2435-2445, Aug. 2008.
<figure><img src="image_23.png" /></figure>

**Lianru Gao** (Senior Member, IEEE) received the B.S. degree in civil engineering from Tsinghua University, Beijing, China, in 2002, and the Ph.D. degree in cartography and geographic information system from the Institute of Remote Sensing Applications, Chinese Academy of Sciences (CAS), Beijing, in 2007.

He is a Professor with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, CAS. He also has been a Visiting Scholar with the University of Extremadura, Cáceres, Spain, in 2014, and Mississippi State University (MSU), Starkville, MS, USA, in 2016. In last 10 years, he was the principal investigator of ten scientific research projects at national and ministerial levels, including projects by the National Natural Science Foundation of China (2010-2012, 2016-2019, 2018-2020), and by the Key Research Program of the CAS (2013-2015) et al. He has authored or coauthored more than 160 peer-reviewed articles, and there are more than 80 journal articles included by Science Citation Index (SCI). He was a coauthor of an academic book *Hyperspectral Image Classification and Target Detection*. He holds 28 National Invention Patents in China. His research focuses on hyperspectral image processing and information extraction. Dr. Gao was awarded the Outstanding Science and Technology Achievement Prize of the CAS in 2016, and was supported by the China National Science Fund for Excellent Young Scholars in 2017, and won the Second Prize of The State Scientific and Technological Progress Award in 2018. He received the recognition of the Best Reviewers of the IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING in 2015, and the Best Reviewers of the IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING in 2017.

<figure><img src="image_24.png" /></figure>

**Haoyang Yu** (S'16-M'19) received the B.S. degree in Information and Computing Science from Northeastern University, Shenyang, China, in 2013, and the Ph.D. degree in cartography and geographic information system from the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences (CAS), Beijing, China, in 2019.

He is currently a Xing Hai Associate Professor with the Center of Hyperspectral Imaging in Remote Sensing (CHIRS), Information Science and Technology College, Dalian Maritime University. His research focuses on models and algorithms for hyperspectral image processing, analysis and applications.

<figure><img src="image_25.png" /></figure>

**Zhicheng Wang** received the B.E. degree from the College of Geosciences and Surveying Engineering, China University of Mining and Technology (Beijing), Beijing, China, in 2018. He is pursuing the M.S. degree in electronics and communication engineering with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences (CAS), Beijing.

His research interests include image processing, machine learning, hyperspectral image nonlinear unmixing, and denoising.

<figure><img src="image_26.png" /></figure>

**Bing Zhang** (Fellow, IEEE) received the B.S. degree in geography from Peking University, Beijing, China, in 1991, and the M.S. and Ph.D. degrees in remote sensing from the Institute of Remote Sensing Applications, Chinese Academy of Sciences (CAS), Beijing, in 1994 and 2003, respectively.

He is a Full Professor and the Deputy Director of the Aerospace Information Research Institute, CAS, where he has been leading lots of key scientific projects in the area of hyperspectral remote sensing for more than 25 years. He has developed five software systems in the image processing and applications. His creative achievements were rewarded ten important prizes from the Chinese government, and special government allowances of the Chinese State Council. He has authored more than 300 publications, including more than 170 journal articles. He has edited six books/contributed book chapters on hyperspectral image processing and subsequent applications. His research interests include the development of Mathematical and Physical models and image processing software for the analysis of hyperspectral remote sensing data in many different areas.

<figure><img src="image_27.png" /></figure>

**Lina Zhuang** (S'15-M'21) received Bachelor's degrees in geographic information system and in economics from South China Normal University, Ghuangzhou, China, in 2012, the M.S. degree in cartography and geography information system from Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing, China, in 2015, and the Ph.D. degree in Electrical and Computer Engineering at the Instituto Superior Tecnico, Universidade de Lisboa, Lisbon, Portugal in 2018.

From 2015 to 2018, she was a Marie Curie Early Stage Researcher of Sparse Representations and Compressed Sensing Training Network (SpaRTaN number 607290) with the Instituto de Telecomunicações. SpaRTaN Initial Training Networks (ITN) is funded under the European Union's Seventh Framework Programme (FP7-PEOPLE-2013-ITN) call and is part of the Marie Curie Actions-ITN funding scheme. From 2019 to 2021, she was a Research Assistant Professor with Hong Kong Baptist University. She is currently a Research Assistant Professor with the University of Hong Kong. Her research interests include hyperspectral image restoration, superresolution, and compressive sensing.

Dr. Zhang has been serving as a Technical Committee Member of IEEE Workshop on Hyperspectral Image and Signal Processing, since 2011, and has been the President of Hyperspectral Remote Sensing Committee of China National Committee of International Society for Digital Earth since 2012, and has been the Standing Director of Chinese Society of Space Research (CSSR) since 2016. He is the Student Paper Competition Committee Member in International Geoscience and Remote Sensing Symposium from 2015 to 2020. He was awarded the National Science Foundation for Distinguished Young Scholars of China in 2013, and the 2016 Outstanding Science and Technology Achievement Prize of the Chinese Academy of Sciences, the highest level of Awards for the CAS scholars. He is serving as an Associate Editor for the IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING.
<figure><img src="image_28.png" /></figure>

**Jocelyn Chanussot** (Fellow, IEEE) received the M.Sc. degree in electrical engineering from the Grenoble Institute of Technology (Grenoble INP), Grenoble, France, in 1995, and the Ph.D. degree from the Université de Savoie, Annecy, France, in 1998.

Since 1999, he has been with Grenoble INP, where he is a Professor of signal and image processing. His research interests include image analysis, hyperspectral remote sensing, data fusion, machine learning and artificial intelligence. He has been a Visiting Scholar with Stanford University, Stanford, CA, USA, KTH Royal Institute of Technology, Stockholm, Sweden, and National University of Singapore, Singapore. Since 2013, he is an Adjunct Professor with the University of Iceland, Reykjavík, Iceland. In 2015–2017, he was a Visiting Professor with the University of California, Los Angeles (UCLA), Los Angeles, CA, USA. He holds the AXA chair in remote sensing and is an Adjunct Professor with the Chinese Academy of Sciences, Aerospace Information Research Institute, Beijing.

Dr. Chanussot is the founding President of IEEE Geoscience and Remote Sensing French chapter (2007–2010) which received the 2010 IEEE Geoscience and Remote Sensing Society Chapter Excellence Award. He has received multiple outstanding paper awards. He was the Vice-President of the IEEE Geoscience and Remote Sensing Society, in charge of meetings and symposia (2017–2019). He was the General Chair of the first IEEE GRSS Workshop on Hyperspectral Image and Signal Processing, Evolution in Remote sensing (WHISPERS). He was the Chair (2009–2011) and Cochair of the GRS Data Fusion Technical Committee (2005–2008). He was a Member of the Machine Learning for Signal Processing Technical Committee of the IEEE Signal Processing Society (2006–2008) and the Program Chair of the IEEE International Workshop on Machine Learning for Signal Processing (2009). He is an Associate Editor for the IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, the IEEE TRANSACTIONS ON IMAGE PROCESSING, and the PROCEEDINGS OF THE IEEE. He was the Editor-in-Chief of the IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING (2011–2015). In 2014, he served as a Guest Editor for the IEEE Signal Processing Magazine. He is a Member of the Institut Universitaire de France (2012–2017) and a Highly Cited Researcher (Clarivate Analytics/Thomson Reuters, 2018–2019).