A Graph Regularized Multilinear Mixing Model for Nonlinear Hyperspectral Unmixing

Abstract

Spectral unmixing of hyperspectral images is an important issue in the fields of remote sensing. Jointly exploring the spectral and spatial information embedded in the data is helpful to enhance the consistency between mixing/unmixing models and real scenarios. This paper proposes a graph regularized nonlinear unmixing method based on the recent multilinear mixing model (MLM). The MLM takes account of all orders of interactions between endmembers, and indicates the pixel-wise nonlinearity with a single probability parameter. By incorporating the Laplacian graph regularizers, the proposed method exploits the underlying manifold structure of the pixels’ spectra, in order to augment the estimations of both abundances and nonlinear probability parameters. Besides the spectrum-based regularizations, the sparsity of abundances is also incorporated for the proposed model. The resulting optimization problem is addressed by using the alternating direction method of multipliers (ADMM), yielding the so-called graph regularized MLM (G-MLM) algorithm. To implement the proposed method on large hypersepectral images in real world, we propose to utilize a superpixel construction approach before unmixing, and then apply G-MLM on each superpixel. The proposed methods achieve superior unmixing performances to state-of-the-art strategies in terms of both abundances and probability parameters, on both synthetic and real datasets.

1. Introduction

Spectral unmixing (SU) has become one of the essential topics in the context of hyperspectral imagery processing. Hyperspectral images are captured by remote sensors over hundreds of channels within a certain wavelength range. Each such image can be viewed as a data cube, with two of its dimensions being the spatial and the third one recording the spectral information, i.e., each pixel corresponds to a spectrum vector [1]. Compared to abundant spectral information, the spatial information contained in the image is quite limited due to the low spatial resolution of sensors. Therefore, it is usually assumed that each observed spectrum is mixed by several pure material signatures, termed endmembers. The aim of SU is to extract the endmembers and to determine their respective proportions, referred to as abundances [2].

The linear mixing model (LMM) is the most widely-used one, which assumes that the spectra mixture occurs on a macroscopic scale. Namely, each photon interacts with only one endmember before being received by the sensor, and the light reflected from different endmembers is then mixed [3]. Therefore, each observed pixel can be represented as a linear combination of the endmembers [4]. To be physically interpretable, the abundance nonnegativity constraint (ANC) and the abundance sum to one constraint (ASC) are often considered in SU algorithms, e.g., the fully constrained least-squares (FCLS) algorithm [5] and the constrained spectral unmixing by splitting and augmented Lagrangian (C-SUnSAL) [6]. A thorough review of the SU based on LMM is provided in [1].

Although LMM-based models can approximate the light reflection and scattering mechanism for most cases, serious nonlinearity effects may exist in many scenes, and more accurate models to describe the mixing/unmixing process of the observed spectra [7] are thus required. To this end, numerous nonlinear unmixing techniques have been developed, including the earlier Hapke’s model [8] for intimate unmixing, the bilinear-based models, and the data-driven, kernel-based unmixing ones [9,10]. The bilinear mixing models essentially assume that each photon interacts with two endmembers successively before reaching the sensor. It is observed that, when the area is covered by soil and plant, bilinear interactions often exist between two materials [11]. Generally, most of the bilinear models are formulated in a general form of adding extra second-order products between endmembers to the LMM. Specific models include the Nascimento model [12], the Fan model [13], the generalized bilinear model (GBM) [14,15], and the polynomial postnonlinear mixing model (PPNMM) [16], and they mainly vary from one to another in the way of imposing constraints on model parameters. In [7], a detailed overview of the bilinear models, as well as their relations and differences, are provided.

Recently, important progress has been made to develop the physics-driven unmixing models accounting for higher order nonlinear effects [17,18,19]. The first model is the so-called multilinear mixing model (MLM) [17]. The MLM extends the bilinear model with the PPNMM formulation to include all degrees of interactions, where a pixel-wise probability parameter is introduced to describe the possibility of light undergoing further interactions. Based on MLM, the authors of [18] proposed an unsupervised nonlinear unmixing method, termed MLMp. The latter improves the former from two aspects: First, the objective function of MLMp is slightly simplified compared to that of MLM. This helps to avoid the under-estimation of the probability parameter in MLM, and to reduce the computational complexity. Second, MLMp is in an unsupervised form, in the sense that the endmembers are estimated jointly with the abundances and probability parameters. In contrast, the MLM is a supervised method, where the endmembers should be extracted in prior by using some endmember estimation algorithms, e.g., vertex component analysis (VCA) [19]. More recently, a novel band-wise nonlinear unmixing algorithm (BNLSU) has been developed in [20]. This model extends the pixel-wise, scalar probability parameter in MLM to the vector form. That is, for a given pixel, its nonlinear intensity at each band is characterized by a probability parameter independently. In particular, a graph-based regularizer is included to enforce similarity between the probability parameters from neighboring bands. In this way, the nonlinear contributions from different bands are expected to be modeled more accurately. Nevertheless, research on MLM model is still very limited.

Jointly exploring the spectral and spatial information contained in the hyperspectral data has shown great effectiveness in boosting the SU performance, mainly by including regularizers favoring different purposes to the objective function of SU models [21,22,23,24,25]. The sparsity constraint, which is originated from sparse unmixing problem [6], consists of expressing the mixing pixels by using only a few endmembers, namely only the abundance elements corresponding to a few endmembers should be non-zero. The sparse solutions are obtained by enforcing $l_0$-norm, or its surrogate $l_1$-norm regularization on the unknown abundances [21,26]. Total-variation (TV) penalties take account of the geographical structure among pixels, by assuming that adjacent pixels tend to have similar abundance vectors for a same set of endmembers. When applied to SU tasks, TV-regularizers usually lead to abundance maps with good spatial consistency [22,27]. By exploiting the intrinsic manifold structure of hyperspectral data, the graph Laplacian regularization has been incorporated within the sparsity-constraint NMF for unsupervised SU [23], and for supervised unmixing [24,28]. Different from TV-inducing regularizers that restrict spatial similarity within some local neighborhood, graph Laplacian regularizer allows one pixel being connected to multiple pixels as long as they are spectrally similar. It is noticed that unmixing methods with graph-based regularizers often suffer from serious computational burden, especially when the hyperspectral image under processing is large. This is because a large graph is built across the entire image, each node stands for a pixel and each edge connects two similar pixels. To alleviate this issue, the authors of [24] firstly applied the spectral clustering algorithm to divide the nodes in the graph into several subgraphs. In [28], the original image is preprocessed using the SLIC-based (simple linear iterative clustering) superpixel construction, by which the adjacent pixels with similar spectral features are divided into small blocks, called superpixels. After these procedures, linear unmixing with graph-based regularization is performed on each subgraph or superpixel with smaller size.

The high-dimensional observed spectra can be taken as being sampled from some sub-manifold [23]. Therefore, it is important to take account of the manifold structure of data in the mixing/unmixing models, in order to augment the parameters estimation. To the best of our knowledge, manifold regularization has not yet been approached for simultaneously regularizing the abundance matrix and the probability parameters for MLM-based methods. In this paper, we propose a novel graph regularized nonlinear SU method based on recent MLM model, and derive the associated algorithm using ADMM. The resulting unmixing algorithm is referred as the graph regularized MLM (G-MLM). The main contributions are as follows.

• By taking advantage of graph Laplacian and sparse unmixing, both the abundance matrix and the nonlinearity parameters are augmented. In particular, the internal manifold structure of data is well-exploited. The intuition is that, if two pixels are spectrally similar with connection established in the graph, their low-rank representations by abundance vectors, as well as their nonlinearity parameters should preserve this consistency.
• We apply the alternating direction method of multipliers (ADMM) [29] to solve the resulting optimization problem. Moreover, to reduce algorithm complexity and improve efficiency, efforts have been made to avoid large graph computation associated to the whole image. To this end, we exploit the superpixel construction method proposed in [28]. By dividing the original image into smaller adjacent superpixels and then performing the unmixing step on each of them, the computational efficiency of the proposed G-MLM algorithm is greatly improved without deteriorating the unmixing performance too much. For the purpose of differentiation, we refer to the graph regularized MLM with superpixel construction by ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$.

The rest of this paper is organized as follows. Section 2 briefly reviews the MLM model and its current variations. In Section 3, the proposed unmixing model is presented, and the optimization problem is formulated. Experimental results on a series of synthetic data and a real image are reported in Section 4. Section 5 concludes the paper with remarks on future works.

2. Related Work

  This section succinctly reviews the existing works related to the MLM unmixing model. Before proceeding, main symbols used in this work are firstly introduced. Given a hyperspectral image, let $\mathit{X}=\left[{\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_T\right]\in {\mathbb{R}}^{L\times T}$ represent the observed data matrix consisting of T pixels over L spectral bands, where ${\mathit{x}}_j\in {\mathbb{R}}^L$ records the spectrum vector of the jth pixel, for $j=1,2,\dots ,T$. Suppose that the hyperspectral data is composed by R endmembers. Let $\mathit{A}=\left[{\mathit{a}}_1,{\mathit{a}}_2,\dots ,{\mathit{a}}_R\right]\in {\mathbb{R}}^{L\times R}$ be the endmember matrix, with ${\mathit{a}}_i\in {\mathbb{R}}^L$ being the spectrum vector of the ith endmember. Let $\mathit{S}=\left[{\mathit{s}}_1,{\mathit{s}}_2,\dots ,{\mathit{s}}_T\right]\in {\mathbb{R}}^{R\times T}$ denote the abundance matrix, where ${\mathit{s}}_j\in {\mathbb{R}}^R$ is the abundance vector of the jth pixel. The entry of abundance matrix, namely $s_{ij}$, denotes the abundance with respect to the ith endmember of the jth pixel. Particularly, we use ${\mathit{y}}_j={\sum }_{i=1}^R{s}_{ij}{\mathit{a}}_i=\mathit{A}{\mathit{s}}_j$ to denote the linear part by LMM of the jth observed pixel ${\mathit{x}}_j$, and $\mathit{Y}=\left[{\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_T\right]\in {\mathbb{R}}^{L\times T}$. Moreover, we use $\mathit{P}=\left[P_1,P_2,\dots ,P_T\right]\in {\mathbb{R}}^{1\times T}$ to denote a row recording the probability parameters for every pixel in the image, where $P_j$ is a scalar indicating the probability of further interactions for the jth pixel.

2.1. MLM

The MLM in [17], a newly-proposed nonlinear mixing model, successfully generalizes the bilinear model to include all degrees of multiple scattering among endmembers. This model is built with a meaningful physical reasoning, upon following basic assumptions:

• The incoming light will interact with at least one material.
• For pixel ${\mathit{x}}_j$, after each interaction with a material, the probability to undergo further interactions is $P_j$, and, accordingly, the probability of escaping the scene and arriving at the sensor becomes $(1-P_j)$.
• The probability of interacting with material i is proportional to its abundance $s_{ij}$.
• The intensity of the light scattered by material i depends on the corresponding material’s albedo ${\mathit{w}}_i\in {\left[0,1\right]}^L$.

As only small differences exist between albedo and reflectance in most practical situations [17], the albedo is substituted by reflectance with $\forall j,{\mathit{w}}_j\approx {\mathit{a}}_j$. Based on these assumptions, the MLM model is formulated as follows:

$$\begin{array}{cc}\hfill {\mathit{x}}_j=& (1-P_j)\sum _{i=1}^R{\mathit{a}}_i{s}_{ij}+(1-P_j)P_j\sum _{i=1}^R\sum _{k=1}^R({\mathit{a}}_i\odot {\mathit{a}}_k)s_{ij}{s}_{kj}\hfill \\ \hfill & +(1-P_j)P_j^2\sum _{i=1}^R\sum _{k=1}^R\sum _{l=1}^R({\mathit{a}}_i\odot {\mathit{a}}_k\odot {\mathit{a}}_l)s_{ij}{s}_{kj}{s}_{lj}+\dots +{\epsilon }_j\hfill \\ \hfill =& (1-P_j){\mathit{y}}_j+P_j{\mathit{y}}_j\odot {\mathit{x}}_j+{\epsilon }_j\hfill \end{array}$$

(1)

where ${\mathit{y}}_j={\displaystyle \sum _{i=1}^R}{s}_{ij}{\mathit{a}}_i=\mathit{A}{\mathit{s}}_j$ denotes the linear part by LMM model, ${\epsilon }_j\in {\mathbb{R}}^L$ is the additive Gaussian noise, and ⊙ represents element-wise product. The resulting optimization problem is given by

$$\begin{array}{cc}\hfill & \underset{{\mathit{s}}_j,P_j}{argmin}{∥{\mathit{x}}_j-\frac{(1-P_j){\mathit{y}}_j}{1-P_j{\mathit{y}}_j}∥}_2^2\hfill \\ \hfill & s.t.{\mathit{s}}_j\ge 0\mathrm{and}{\mathbf{1}}_R^{\top }{\mathit{s}}_j=1\hfill \\ \hfill & P_j\le 1\hfill \end{array}$$

(2)

where ${\mathbf{1}}_R\in {\mathbb{R}}^R$ represents the vector of ones. The sequential quadratic programming is applied to solve above optimization problem. Note that the MLM is a supervised unmixing model in the sense that the endmember matrix $\mathit{A}$ is extracted by VCA [19] in prior. The pixel-wise parameter $P_j$ indicates the probability for further interactions, where a particular case with $P_j=0$ leads to the LMM model. Of particular note is that the model in Equation (1) is also well defined for $P_j<0$.

2.2. Unsupervised MLM (MLMp)

In [18], an unsupervised unmixing method is developed based on MLM. The optimization problem in the so-called MLMp is formulated by

$$\begin{array}{cc}\hfill & \underset{\mathit{A},{\{{\mathit{s}}_j,P_j\}}_{j=1}^T}{argmin}\sum _{j=1}^T{∥(1-P_j){\mathit{y}}_j+P_j{\mathit{y}}_j\odot {\mathit{x}}_j-{\mathit{x}}_j∥}_2^2\hfill \\ \hfill & s.t.{\mathit{s}}_j\ge 0\mathrm{and}{\mathbf{1}}_R^{\top }{\mathit{s}}_j=1\hfill \\ \hfill & P_j\le 1\hfill \\ \hfill & {\mathit{y}}_j=\mathit{A}{\mathit{s}}_j\hfill \\ \hfill & 0\le \mathit{A}\le 1\hfill \end{array}$$

(3)

Note that $0\le \mathit{A}\le 1$ is element-wise. Compared to MLM, MLMp has a simplified objective function without the denominator. This helps to alleviate the under-estimation issue of probability parameters in MLM to some extent, and also to decrease the difficulty in solving the problem in Equation (3) by using a block coordinate descent (BCD) strategy. Another merit is that MLMp is an unsupervised method, which estimates the endmember matrix jointly with the abundances and the probability parameters.

2.3. BNLSU

More recently, the authors of [20] extended MLM for the supervised band-wise nonlinear unmixing (BNLSU). Taking account of the wavelength dependent multiple scatterings, every pixel is supposed to have a specific probability vector recording the band-wise nonlinearities. Considering all the pixels of an image, a probability matrix $\overline{\mathit{P}}=\left[{\overline{\mathit{p}}}_1,{\overline{\mathit{p}}}_2,\dots ,{\overline{\mathit{p}}}_T\right]\in {\mathbb{R}}^{L\times T}$ is introduced, where each vector ${\overline{\mathit{p}}}_j\in {\mathbb{R}}^L$, for $j=1,2,\dots ,T$ describes the probability vector for the jth pixel. In addition, both the sparsity constraint on abundances and the smoothness constraint on probability parameters are incorporated to the BNLSU model, yielding

$$\begin{array}{cc}\hfill \underset{\mathit{S},\overline{\mathit{P}}}{argmin}& \frac{1}{2}\sum _{j=1}^T{∥({\mathbf{1}}_n-{\overline{\mathit{p}}}_j)\odot {\mathit{y}}_j+{\overline{\mathit{p}}}_j{\mathit{y}}_j\odot {\mathit{x}}_j-{\mathit{x}}_j∥}_F^2+{\lambda }_1{\parallel \mathit{S}\parallel }_{1,1}+\frac{{\lambda }_2}{2}\mathrm{Tr}\left({\overline{\mathit{PLP}}}^{\top }\right)\hfill \\ \hfill s.t.& {\mathit{s}}_j\ge 0\mathrm{and}{\mathbf{1}}_R^{\top }{\mathit{s}}_j=1\hfill \\ \hfill & {\overline{\mathit{p}}}_j\le \mathbf{1}\hfill \end{array}$$

(4)

where ${\lambda }_1$ and ${\lambda }_2$ are regularization parameters balancing the influences of constraints, $\mathrm{Tr}(·)$ is the trace of matrix, the inequality ${\overline{\mathit{p}}}_j\le \mathbf{1}$ is element-wise, $\overline{\mathit{L}}=\overline{\mathit{D}}-\overline{\mathit{W}}$ is the graph Laplacian matrix and $\overline{\mathit{D}}$ is a diagonal matrix with the entries $\overline{{\mathit{D}}_{ii}}={\displaystyle \sum _{j=1}^L}\overline{{\mathit{W}}_{ij}}$. Here, the weight matrix $\overline{\mathit{W}}\in {\mathbb{R}}^{L\times L}$ is constructed by measuring the similarity between different pairs of bands across the image. In this work, the Laplacian graph regularizer is adopted to promote smoothness on the estimated probability matrix $\overline{\mathit{P}}$. The optimization problem in Equation (4) is addressed by ADMM. Owing to an increasing number of probability parameters, the BNLSU model shows great effectiveness in reducing the reconstruction error at each band. However, BNLSU is very sensitive to noise due to the calculation of a large number of probability parameters.

3. Proposed Graph Regularized Multilinear Mixing Model

  In this section, we propose a novel graph regularized multilinear mixing model for nonlinear unmixing problem. We first present the constraints utilized in the proposed model, and formulate the optimization problem. Next, we derive the algorithm for solving the resulting optimization problem by applying the ADMM.

3.1. Constraints and Problem Formulation

The abundance sparsity is widely considered in the context of hyperspectral unmixing problem [30,31]. For these cases, the number of endmembers participating in the mixing process is small, leading to the sparse abundance matrix $\mathit{S}$ containing many values of zero. Traditionally, $l_0$ norm can be used to represent the number of nonzero entries in a matrix. Since the resulting optimizations are hard, $l_1$ norm is usually considered instead of $l_0$ for sparsity regularization, with 

$${\parallel \mathit{S}\parallel }_{1,1}=\sum _{j=1}^T{\parallel {\mathit{s}}_j\parallel }_1$$

(5)

Besides the commonly-used sparsity constraint on abundance, the proposed model mainly exploits the underlying manifold structure of the observed spectra, by jointly including the manifold-based regularizations not only on abundance vectors but also on nonlinear probability parameters. Inspired by the work in [23,24], we adopt the Laplacian graph regularizer to promote similarities between similar pixels. To this end, the input image $\mathit{X}\in {\mathbb{R}}^{L\times T}$ is firstly mapped to a graph $\mathit{G}$, where each node of the graph represents a L-dimensional pixel spectrum ${\mathit{x}}_j$, for $j=1,2,\dots ,T$. The affinity matrix $\mathit{W}$ corresponding to the graph should satisfy following requirement: Its element ${\mathit{W}}_{ij}$ reflects the similarity level between pixels ${\mathit{x}}_i$ and ${\mathit{x}}_j$. Specifically, the value of ${\mathit{W}}_{ij}$ is positively correlated to the similarity level. The value decreases to zero for dissimilar pixel pairs. There are different manners to define $\mathit{W}$, e.g., by using heat kernel as in [23,24]. In this paper, we define the affinity matrix $\mathit{W}$ of the graph simply by

$${\mathit{W}}_{ij}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\parallel {\mathit{x}}_i-{\mathit{x}}_j{\parallel }^2<d_{min}^2,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(6)

where $d_{min}^2$ denotes a pre-defined threshold for the squared distance between pixels [24]. It is natural to suppose that, if two pixels ${\mathit{x}}_i$ and ${\mathit{x}}_j$ are spectrally similar, their low-rank representations with respect to abundance vectors ${\mathit{s}}_i$ and ${\mathit{s}}_j$, and with respect to nonlinear probability parameters $P_i$ and $P_j$, should both retain this consistency in respective new spaces. Similar to in [23], by using the affinity matrix $\mathit{W}$, the above assumption on abundances consistency between similar pixels can be modeled by the following regularization

$$\begin{array}{cc}\hfill \frac{1}{2}\sum _{i,j=1}^T{\parallel {\mathit{s}}_i-{\mathit{s}}_j\parallel }^2{\mathit{W}}_{ij}& =\sum _{i=1}^T{\mathit{s}}_i^{\top }{\mathit{s}}_i{\mathit{D}}_{ii}-\sum _{i,j=1}^T{\mathit{s}}_i^{\top }{\mathit{s}}_j{\mathit{W}}_{ij}\hfill \\ \hfill & =\mathrm{Tr}\left(\mathit{S}\mathit{D}{\mathit{S}}^{\top }\right)-\mathrm{Tr}\left(\mathit{S}\mathit{W}{\mathit{S}}^{\top }\right)\hfill \\ \hfill & =\mathrm{Tr}\left(\mathit{S}\mathit{L}{\mathit{S}}^{\top }\right)\hfill \end{array}$$

(7)

where $\mathit{D}$ is a diagonal matrix with the entries ${\mathit{D}}_{ii}={\displaystyle \sum _{j=1}^T}{\mathit{W}}_{ij}$, and $\mathit{L}$ is the graph Laplacian matrix deduced by $\mathit{L}=\mathit{D}-\mathit{W}$.

Analogously, based on the assumption that spectrally similar pixels tend to have nonlinear probability parameters close to each other, the regularization for nonlinear probability parameters is expressed by

$$\frac{1}{2}\sum _{i,j=1}^T{\parallel P_i-P_j\parallel }^2{\mathit{W}}_{ij}=\mathrm{Tr}\left(\mathit{P}\mathit{L}{\mathit{P}}^{\top }\right)$$

(8)

We propose a supervised unmixing method based on the MLM, where the endmembers are supposed to be extracted in prior using some endmember extraction strategy, e.g., VCA [19] and N-Findr [32]. The objective function in Equation (3) in MLMp [18] is taken into consideration in this paper. After incorporating the constraints in Equations (5)–(8) discussed previously, the proposed unmixing model is defined by considering the following optimization problem

$$\begin{array}{cc}\hfill \underset{\mathit{S},\mathit{P}}{argmin}& \frac{1}{2}\sum _{j=1}^T{∥(1-P_j){\mathit{y}}_j+P_j{\mathit{y}}_j\odot {\mathit{x}}_j-{\mathit{x}}_j∥}_F^2\hfill \\ \hfill & +{\lambda }_1{\parallel \mathit{S}\parallel }_{1,1}+\frac{{\lambda }_2}{2}\mathrm{Tr}\left(\mathit{S}\mathit{L}{\mathit{S}}^{\top }\right)+\frac{{\lambda }_3}{2}\mathrm{Tr}\left(\mathit{P}\mathit{L}{\mathit{P}}^{\top }\right)\hfill \\ \hfill s.t.& {\mathit{s}}_j\ge 0,{\mathbf{1}}_R^{\top }{\mathit{s}}_j=1\hfill \\ \hfill & P_j\le 1\hfill \end{array}$$

(9)

where ${\lambda }_1$, ${\lambda }_2$ and ${\lambda }_3$ are hyper-parameters which balance the importances of different regularization terms.

3.2. Nonlinear Unmixing Using ADMM

We utilize the well-known ADMM method [29] to address the optimization problem in Equation (9). The ADMM aims to decompose a hard optimization problem into a sequence of simpler subproblems, and conquers them in an alternating manner. To this end, we introduce new variables $\mathit{G}$ and $\mathit{H}$ for $\mathit{S}$ and $\mathit{P}$, respectively, and reformulate the optimization problem in Equation (9) with

$$\begin{array}{cc}\hfill \underset{\mathit{S},\mathit{P},\mathit{G},\mathit{H}}{argmin}& \frac{1}{2}\sum _{j=1}^T{∥(1-P_j){\mathit{y}}_j+P_j{\mathit{y}}_j\odot {\mathit{x}}_j-{\mathit{x}}_j∥}_F^2+{\lambda }_1{\parallel \mathit{G}\parallel }_{1,1}+\frac{{\lambda }_2}{2}\mathrm{Tr}\left(\mathit{S}\mathit{L}{\mathit{S}}^{\top }\right)+\frac{{\lambda }_3}{2}\mathrm{Tr}\left(\mathit{H}\mathit{L}{\mathit{H}}^{\top }\right)\hfill \\ \hfill & +{\iota }_{\left\{{\mathbf{1}}_T^{\top }\right\}}\left({\mathbf{1}}_R^{\top }\mathit{S}\right)+{\iota }_{{\mathbb{R}}_{+}^{R\times T}}\left(\mathit{G}\right)+{\iota }_{\left\{\mathit{H}\right|\mathit{H}\le {\mathbf{1}}_T^{\top }\}}\left(\mathit{H}\right)\hfill \\ \hfill s.t.& \mathit{G}=\mathit{S},\mathit{H}=\mathit{P}\hfill \end{array}$$

(10)

where ${\iota }_{\left\{{\mathbf{1}}_T^{\top }\right\}}\left({\mathbf{1}}_R^{\top }\mathit{S}\right)$, ${\iota }_{{\mathbb{R}}_{+}^{R\times T}}\left(\mathit{G}\right)$ and ${\iota }_{\left\{\mathit{H}\right|\mathit{H}\le {\mathbf{1}}_T^{\top }\}}\left(\mathit{H}\right)$ are three indicator functions to project ${\mathbf{1}}_R^{\top }\mathit{S}$, $\mathit{G}$, and $\mathit{H}$ onto the sets $\left\{{\mathbf{1}}_R^{\top }\mathit{S}\right|{\mathbf{1}}_R^{\top }\mathit{S}={\mathbf{1}}_T^{\top }\}$, $\left\{\mathit{G}\right|\mathit{G}\ge 0\}$, and $\left\{\mathit{H}\right|\mathit{H}\le {\mathbf{1}}_T^{\top }\}$. Specifically, the ASC and ANC constraints are imposed on $\mathit{S}$ and $\mathit{G}$, respectively, while the constraint on nonlinearity vector is imposed on $\mathit{H}$.

The augmented Lagrangian of Equation (10) can be easily derived by

$$\begin{array}{cc}\hfill \underset{\mathit{S},\mathit{P},\mathit{G},\mathit{H}}{argmin}& \frac{1}{2}\sum _{j=1}^T{∥(1-P_j){\mathit{y}}_j+P_j{\mathit{y}}_j\odot {\mathit{x}}_j-{\mathit{x}}_j∥}_F^2+{\lambda }_1{\parallel \mathit{G}\parallel }_{1,1}+\frac{{\lambda }_2}{2}\mathrm{Tr}\left(\mathit{S}\mathit{L}{\mathit{S}}^{\top }\right)+\frac{{\lambda }_3}{2}\mathrm{Tr}\left(\mathit{H}\mathit{L}{\mathit{H}}^{\top }\right)\hfill \\ \hfill & +{\iota }_{\left\{{\mathbf{1}}_T^{\top }\right\}}\left({\mathbf{1}}_R^{\top }\mathit{S}\right)+{\iota }_{{\mathbb{R}}_{+}^{R\times T}}\left(\mathit{G}\right)+{\iota }_{\left\{\mathit{H}\right|\mathit{H}\le {\mathbf{1}}_T^{\top }\}}\left(\mathit{H}\right)+\frac{\rho }{2}({∥\mathit{S}-\mathit{G}+{\mathit{M}}_1∥}_F^2+{∥\mathit{P}-\mathit{H}+{\mathit{M}}_2∥}_F^2)\hfill \end{array}$$

(11)

where $\rho$ is the penalty parameter which is usually set to be a small positive value, and ${\mathit{M}}_1$ and ${\mathit{M}}_2$ are the scaled dual variables, which are of the same size as $\mathit{S}$ and $\mathit{P}$, respectively.

At each iteration, the ADMM algorithm minimizes the augmented Lagrangian in Equation (11) iteratively, by alternating the minimization over each of the variables. Namely, every variable is alternately updated while keeping the other variables fixed to their latest values. For each iteration, we start with the minimization with respect to $\mathit{S}$ and $\mathit{G}$, and then with respect to $\mathit{P}$ and $\mathit{H}$. Finally, we update the scaled dual variables ${\mathit{M}}_1$ and ${\mathit{M}}_2$. During the optimization, the primal residual

$${\mathrm{res}}_{\mathrm{p}}={∥\begin{array}{c}\mathit{S}-\mathit{G}\\ \mathit{P}-\mathit{H}\end{array}∥}_F$$

(12)

and the dual residual

$${\mathrm{res}}_{\mathrm{d}}={∥\begin{array}{c}{\mathit{G}}_0-\mathit{G}\\ {\mathit{H}}_0-\mathit{H}\end{array}∥}_F$$

(13)

are examined over iterations, in order to check whether the stopping-condition is attained, where ${\mathit{G}}_0$ and ${\mathit{H}}_0$ represent the values in previous iteration.

3.2.1. Update S and G

At the $(t+1)$th iteration, we update $\mathit{S}$ to ${\mathit{S}}^{(t+1)}$. By discarding the terms irrelevant with $\mathit{S}$ in Equation (11), the reduced optimization subproblem becomes

$${\mathit{S}}^{(t+1)}=\underset{\mathit{S}}{argmin}\frac{1}{2}\sum _{j=1}^T{∥(1-P_j){\mathit{y}}_j+P_j{\mathit{y}}_j\odot {\mathit{x}}_j-{\mathit{x}}_j∥}_F^2+{\iota }_{\left\{1\right\}}\left({\mathbf{1}}_R^{\top }{\mathit{s}}_j\right)+\frac{{\lambda }_2}{2}\mathrm{Tr}\left(\mathit{S}\mathit{L}{\mathit{S}}^{\top }\right)+\frac{\rho }{2}{∥\mathit{S}-{\mathit{G}}^{\left(t\right)}+{\mathit{M}}_1^{\left(t\right)}∥}_F^2$$

(14)

Let ${\overline{{\mathit{A}}_j}}^{\left(t\right)}=\mathit{A}\odot \left((1-P_j^{\left(t\right)}){\mathbf{1}}_{L\times R}+P_j^{\left(t\right)}{\mathit{x}}_j{\mathbf{1}}_R^{\top }\right)$ be the pixel-wise matrix depending on $P_j^{\left(t\right)}$, for $j=1,2,\dots T$. Then, the problem in Equation (14) becomes

$${\mathit{S}}^{(t+1)}=\underset{\mathit{S}}{argmin}\frac{1}{2}\sum _{j=1}^T{∥{\overline{{\mathit{A}}_j}}^{\left(t\right)}{\mathit{s}}_j-{\mathit{x}}_j∥}_F^2+{\iota }_{\left\{1\right\}}\left({\mathbf{1}}_R^{\top }{\mathit{s}}_j\right)+\frac{{\lambda }_2}{2}\mathrm{Tr}\left(\mathit{S}\mathit{L}{\mathit{S}}^{\top }\right)+\frac{\rho }{2}{∥\mathit{S}-{\mathit{G}}^{\left(t\right)}+{\mathit{M}}_1^{\left(t\right)}∥}_F^2,$$

(15)

or equivalently in the following vector-wise form

$${\mathit{s}}_j^{(t+1)}=\underset{{\mathit{s}}_j}{argmin}\frac{1}{2}{∥{\overline{\mathit{A}}}_j^{\left(t\right)}{\mathit{s}}_j-{\mathit{x}}_j∥}_F^2+\frac{1}{2}\sum _{l=1}^T{\parallel {\mathit{s}}_l-{\mathit{s}}_j\parallel }^2{\mathit{W}}_{ij}+{\iota }_{\left\{1\right\}}\left(1^{\top }{\mathit{s}}_j\right)+\frac{\rho }{2}{∥{\mathit{s}}_j-{\mathit{g}}_j^{\left(t\right)}+{\left({\mathit{m}}_1\right)}_j^{\left(t\right)}∥}_F^2,$$

(16)

for $j=1,2,\dots ,T$. Similar to in [6], the solution of Equation (16) is derived by solving a quadratic problem with equality constraints, given by

$${\mathit{s}}_j^{(t+1)}={\mathit{B}}^{-1}\mathit{w}-\mathit{C}({\mathbf{1}}_R^{\top }{\mathit{B}}^{-1}\mathit{w}-{\mathbf{1}}_R),$$

(17)

with

$$\begin{array}{ccc}\hfill \mathit{B}& =& {\left({\overline{\mathit{A}}}_j^{\left(t\right)}\right)}^{\top }{\overline{\mathit{A}}}_j^{\left(t\right)}+(\rho +\frac{{\lambda }_2}{2}\sum _{l=1}^T{\mathit{W}}_{lj}){\mathit{I}}_{R\times R}\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill \mathit{C}& =& {\mathit{B}}^{-1}{\mathbf{1}}_R{\left({\mathbf{1}}_R^T{\mathit{B}}^{-1}{\mathbf{1}}_R\right)}^{-1}\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill \mathit{w}& =& {\left({\overline{\mathit{A}}}_j^{\left(t\right)}\right)}^{\top }{\mathit{x}}_j+\rho \left({\mathit{g}}_j^{\left(t\right)}-{\left({\mathit{m}}_1\right)}_j^{\left(t\right)}\right)+\frac{{\lambda }_2}{2}\sum _{l=1}^T{\mathit{W}}_{lj}{\mathit{s}}_l\hfill \end{array}$$

(20)

where ${\mathit{g}}_j^{\left(t\right)}$ and ${\left({\mathit{m}}_1\right)}_j^{\left(t\right)}$ denote the columns of $\mathit{G}$ and ${\mathit{M}}_1$, respectively, and ${\mathit{I}}_{R\times R}$ is the identity matrix.

Regarding the optimization in terms of $\mathit{G}$, the corresponding reduced subproblem is

$${\mathit{G}}^{(t+1)}=\underset{\mathit{G}}{argmin}{\lambda }_1{\parallel \mathit{G}\parallel }_{1,1}+{\iota }_{R_{+}^{R\times T}}\left(\mathit{G}\right)+\frac{\rho }{2}{∥{\mathit{S}}^{(t+1)}-\mathit{G}+{\mathit{M}}_1^{\left(t\right)}∥}_F^2$$

(21)

If we ignore the ANC constraint enforced by ${\iota }_{R_{+}^{R\times T}}\left(\mathit{G}\right)$, the solution of Equation (21) would be

$${\mathit{G}}^{(t+1)}={\mathcal{S}}_{{\lambda }_1/\rho }({\mathit{S}}^{(t+1)}+{\mathit{M}}_1^{\left(t\right)})$$

(22)

where $S_{{\lambda }_1/\rho }$ is the soft thresholding operator [29] defined by

$${\mathcal{S}}_b\left(\zeta \right)=\left\{\begin{array}{ccc}\hfill \zeta -b& & \zeta >b\hfill \\ \hfill 0& & \left|\zeta \right|<b\hfill \\ \hfill \zeta +b& & \zeta <-b\hfill \end{array}\right.$$

(23)

To further impose the ANC constraint, it is straightforward to project the result in Equation (21) onto the first quadrant according to [6], yielding

$${\mathit{G}}^{(t+1)}=max\{{\mathbf{0}}_{R\times T},{\mathit{G}}^{(t+1)}\}$$

(24)

where the maximum function is applied in an element-wise manner.

3.2.2. Update P and H

By discarding the terms independent of $\mathit{P}$, the reduced optimization problem of $\mathit{P}$ is expressed by

$${\mathit{P}}^{(t+1)}=\underset{\mathit{P}}{argmin}\frac{1}{2}\sum _{j=1}^T{∥(1-P_j){\mathit{y}}_j+P_j{\mathit{y}}_j\odot {\mathit{x}}_j-{\mathit{x}}_j∥}_F^2+\frac{\rho }{2}{∥\mathit{P}-{\mathit{H}}^{\left(t\right)}+{\mathit{M}}_2^{\left(t\right)}∥}_F^2$$

(25)

The update rule computes

$${\mathit{P}}^{(t+1)}=\left({\mathbf{1}}_L^{\top }\left(({\mathit{Y}}^{(t+1)}-{\mathit{Y}}^{(t+1)}\odot \mathit{X})\odot ({\mathit{Y}}^{(t+1)}-\mathit{X})\right)+\rho \left({\mathit{H}}^{\left(t\right)}-{\mathit{M}}_2^{\left(t\right)}\right)\right)./\left({\mathbf{1}}_L^{\top }{({\mathit{Y}}^{(t+1)}-{\mathit{Y}}^{(t+1)}\odot \mathit{X})}^2+\rho {\mathbf{1}}_T^{\top }\right),$$

(26)

where ${\mathit{Y}}^{(t+1)}=\mathit{A}{\mathit{S}}^{(t+1)}$ is a matrix recording LMM terms as previously defined, and both $./$ and ${(·)}^2$ indicate that the operations are element-wise.

In a similar way, the subproblem with respect to $\mathit{H}$ is given by

$${\mathit{H}}^{(t+1)}=\underset{\mathit{H}}{argmin}\frac{{\lambda }_3}{2}\mathrm{Tr}\left(\mathit{H}\mathit{L}{\mathit{H}}^{\top }\right)+{\iota }_{\left\{\mathit{H}\right|\mathit{H}\le {\mathbf{1}}_T^{\top }\}}\left(\mathit{H}\right)+\frac{\rho }{2}{∥{\mathit{P}}^{(t+1)}-\mathit{H}+{\mathit{M}}_2^{\left(t\right)}∥}_F^2$$

(27)

Without the boundary constraint ${\iota }_{\left\{\mathit{H}\right|\mathit{H}\le {\mathbf{1}}_T^{\top }\}}\left(\mathit{H}\right)$ on nonlinearity parameters, the solution of Equation (27) would be

$${\mathit{H}}^{(t+1)}=\rho ({\mathit{P}}^{(t+1)}+{\mathit{M}}_2^{\left(t\right)}){({\lambda }_3\mathit{L}+\rho {\mathbf{1}}_{T\times T})}^{-1}$$

(28)

By adding the boundary constraint, the update rule of ${\mathit{H}}^{(t+1)}$ becomes

$${\mathit{H}}^{(t+1)}=min\{{\mathbf{1}}_T^{\top },{\mathit{H}}^{(t+1)}\}$$

(29)

where the minimum function is element-wise.

3.2.3. Update $M_1$ and $M_2$

The last step of ADMM consists of updating the scaled dual variables ${\mathit{M}}_1$ and ${\mathit{M}}_2$ in the following manner

$${\mathit{M}}_1^{(t+1)}={\mathit{M}}_1^{\left(t\right)}+{\mathit{S}}^{(t+1)}-{\mathit{G}}^{(t+1)}$$

(30)

$${\mathit{M}}_2^{(t+1)}={\mathit{M}}_2^{\left(t\right)}+{\mathit{P}}^{(t+1)}-{\mathit{H}}^{(t+1)}$$

(31)

where ${\mathit{M}}_1$ and ${\mathit{M}}_2$ measure the sum of the primal residuals, reflecting the compliance levels of $\mathit{S}$ and $\mathit{P}$ with their constraints, respectively.

According to Boyd et al. [29], a good initialization is often beneficial to convergence of ADMM. In this paper, we use the fully constrained least squares (FCLS) [5] to initialize the abundance matrix ${\mathit{S}}^{\left(0\right)}$, such that the result satisfies both ANC and ASC. On this basis, the nonlinearity parameters vector ${\mathit{P}}^{\left(0\right)}$ is initialized according to Equation (26), where the boundary constraint ${\mathit{P}}^{\left(0\right)}\le 1$ is satisfied. Let ${\mathit{G}}^{\left(0\right)}$ and ${\mathit{H}}^{\left(0\right)}$ be equal to ${\mathit{S}}^{\left(0\right)}$ and ${\mathit{P}}^{\left(0\right)}$, and ${\mathit{M}}_1^{\left(0\right)}$ and ${\mathit{M}}_2^{\left(0\right)}$ be initialized by zero matrix and zero vector, respectively. The stopping condition of Algorithm 1 for G-MLM is two-fold: (1) The primal residual in Equation (12) and dual residual in Equation (13) are smaller than the pre-defined tolerance values, namely ${\mathrm{res}}_{\mathrm{p}}\le \mathrm{tol}1$ and ${\mathrm{res}}_{\mathrm{d}}\le \mathrm{tol}2$. In the experiments, we set $\mathrm{tol}1=\mathrm{tol}2=\sqrt{\mathrm{T}\mathrm{R}}\times {10}^{-4}$, following [6]. (2) The maximum number of iterations does not exceed the preset value ${\mathrm{Iter}}_{\max }$.

We analyze the computational complexity of the proposed G-MLM in Algorithm 1, by specifying the updating complexity for each variable, as shown in  

Table 1. The complexity of the proposed G-MLM for updating each variable per iteration.

Variable to Update	S	G	P	H	${\mathit{M}}_1$	${\mathit{M}}_2$
Complexity per Iteration	$\mathcal{O}\left(T^{2}LR\right)$	$\mathcal{O}\left(TR\right)$	$\mathcal{O}\left(TL\right)$	$\mathcal{O}\left(T^2\right)$	$\mathcal{O}\left(TR\right)$	$\mathcal{O}\left(T\right)$

. For each iteration, it is the $\mathit{S}$-update and $\mathit{H}$-update that are the most costly ones, as their complexities are $\mathcal{O}\left(T^{2}LR\right)$ and $\mathcal{O}\left(T^2\right)$, with T standing for the total number of pixels. To alleviate the computational burden when tackling large image, we exploit a preprocessing approach before applying G-MLM.

Algorithm 1 ADMM for the proposed graph regularized MLM (G-MLM)

Input:$\mathit{X}\in {\mathbb{R}}^{L\times T}$: hyperspectral data; $\mathit{A}$: endmember matrix  1:   Output:$\mathit{S}\in {\mathbb{R}}^{R\times T}$: abundance matrix; $\mathit{P}\in {\mathbb{R}}^{1\times T}$: probability parameter vector  2:    3: set$t=0$, $\rho >0$, ${\lambda }_1>0$, ${\lambda }_2>0$ and ${\lambda }_3>0$.  4:    5: initialize${\mathit{S}}^{\left(0\right)}$, ${\mathit{P}}^{\left(0\right)}$, ${\mathit{G}}^{\left(0\right)}$, ${\mathit{H}}^{\left(0\right)}$, ${\mathit{M}}_1^{\left(0\right)}$, and ${\mathit{M}}_2^{\left(0\right)}$.  6:    7: while stopping conditions are not met, do  8:    9:     update ${\mathit{S}}^{(t+1)}$ and ${\mathit{G}}^{(t+1)}$ with Equations (17) and (24), respectively. 10:   11:     update ${\mathit{P}}^{(t+1)}$ and ${\mathit{H}}^{(t+1)}$ with Equations (26) and (29), respectively. 12:   13:     update ${\mathit{M}}_1^{(t+1)}$ and ${\mathit{M}}_2^{(t+1)}$ with Equations (30) and (31), respectively. 14:   15:     $t=t+1$ 16:   17: end while

• The first strategy is to use the spectral clustering algorithm proposed in [33], following Ammanouil et al. [24]. This method divides the nodes of the original graph into k subgraphs, and the affinity matrix of each subgraph consists of a subset of the original affinity matrix. Then, the proposed G-MLM is performed on every subgraph, which contains a number of pixels far less than T. We refer the G-MLM algorithm combined with spectral clustering method [24] by ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$. To balance the computational complexity and the global knowledge represented by the graph, the number of subgraphs k should be set conservatively, as indicated in [24]. As shown in the Experimental Section 4, although ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ can reduce the computational cost to some extent, it still requires large memory to establish the connections between all pixel pairs across the image.
• The second scheme is to apply a SLIC-based (simple linear iterative clustering) superpixel construction method [28], and to perform G-MLM on each of the superpixels. The superpixel construction in [28] consists of dividing the original image into m small and adjacent superpixels of irregular shapes, according to the spectral similarity of pixels. After that, the G-MLM unmixing method is directly performed on each superpixel that are of much smaller size compared to the original image. In this paper, the resulting method is termed as ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$. As demonstrated next, ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ is effective in reducing the computational time when addressing large image, without deteriorating the unmixing performance of G-MLM too much.

4. Experimental Results

In this section, the proposed G-MLM method, as well as ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ and ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$, is compared with six state-of-the-art unmixing methods, on both synthetic and real hyperspectral datasets. The compared unmixing methods include both the linear and nonlinear ones, namely the C-SUnSAL [6] dedicated to the linear mixing model, the so-called GBM [14], PPNMM [16] elaborated for the bilinear mixing model, and the aforementioned MLM [17], MLMp [18], and BNLSU [20] based on the multilinear mixing model. All algorithms were executed using MATLAB R2016a on a computer with Intel Core i5-8250 CPU at 1.60 GHz and 8GB RAM.

We adopedt following metrics to evaluate the unmixing performance of all the comparing methods. For synthetic datasets, as the ground truth abundance matrix is available, the unmixing result was evaluated by the the root-mean-square error (RMSE) of abundance, defined by

$$\mathrm{RMSE}=\sqrt{\frac{1}{TR}\sum _{j=1}^T{\parallel {\mathit{s}}_j-\overline{{\mathit{s}}_j}\parallel }_2^2},$$

(32)

where ${\mathit{s}}_j$ and $\overline{{\mathit{s}}_j}$ stand for the actual and estimated abundance vector for the jth pixel, respectively. Concerning the real hyperspecral data, as the actual abundance matrix is unknown, the reconstruction error (RE) of pixels was considered for assessing unmixing performance, which is given by

$$\mathrm{RE}=\sqrt{\frac{1}{TL}\sum _{j=1}^T{\parallel {\mathit{x}}_j-\overline{{\mathit{x}}_j}\parallel }_2^2},$$

(33)

where ${\mathit{x}}_j$ and $\overline{{\mathit{x}}_j}$ represent the observed spectrum and reconstructed one for the jth pixel, respectively.

4.1. Experiments on Synthetic Data Cubes

Experiments were firstly conducted on two sets of synthetic data cubes, namely DC1 and DC2 [21,22,34], that were generated by MLM [17] and PPNMM [16], respectively. For either set, there are three images, and each has a size of $75\times 75$ pixels. There are five endmembers for data generation, i.e., $[{\mathit{a}}_1,{\mathit{a}}_2,{\mathit{a}}_3,{\mathit{a}}_4,{\mathit{a}}_5]$, that were randomly selected from the USGS library [35], each one being measured over 224 bands, as shown in Figure 1. The generated abundance maps corresponding to the five endmembers are illustrated in the first line of 5.

There are 25 regions of square in the image, each of the size $5\times 5$ pixels. Each square is homogeneous in terms of fractional abundances, namely the pixels in a same square have the same abundance vectors. The five squares of the first row are pure regions, each composed by only one endmember. For example, the abundance vectors for the square located at the first row and first column are $\mathit{s}={\left[10000\right]}^{\top }$. The five squares of the second row are mixed by two of the endmembers with the same proportions. Taking the square located at the second row and first column as an example, the corresponding abundance vectors are $\mathit{s}={\left[0.50000.5\right]}^{\top }$. Furthermore, the five squares in the last row of the image is a mixture of five endmembers with equal fractions, therefore they are identical. Apart from the 25 squares, the rest pixels in the image are all considered as background, which were generated using the following abundance vector ${\left[0.11490.07410.20030.20550.4051\right]}^{\top }$. Generally, every synthetic image was generated using the aforementioned endmembers and abundance maps, according to different mixture models. The mixing models, as well as the associated parameters for generating DC1 and DC2, are detailed as follows:

DC1: 

The pixels in DC1 images are generated by MLM in Equation (1). Specifically, the nonlinear probability parameter P is distributed as follows: the p values for the background pixels are set to zero by assuming that the backgrounds are linearly mixed. Concerning the remaining 25 regions in the square, each one is assigned by an identical p value, which is randomly chosen from a half-normal distribution with $\sigma =0.3$. The values of p larger than one are set to zero as in [17]. It is also noteworthy that, since the five squares in the last row are identical, the same p value is assigned to all of them.

DC2: 

The synthetic pixels in DC2 images are generated based on the PPNMM [16], which is given by

$${\mathit{x}}_j=\sum _{i=1}^R{s}_{ij}{\mathit{a}}_i+b\left(\sum _{i=1}^R{s}_{ij}{\mathit{a}}_i\right)\odot \left(\sum _{i=1}^R{s}_{ij}{\mathit{a}}_i\right)+{\epsilon }_j,$$

(34)

where the nonlinear parameter b is generated similarly as in DC1. For the assemblage of background pixels and the remaining 25 squares, each is assigned by the same b randomly drawn from a uniform distribution within the range $[-0.3,0.3]$, according to Altmann et al. [16].

Finally, to test the robustness of the proposed G-MLM against noise, the Gaussian noise of varying signal-to-noise ratio (SNR) is added, with $\mathrm{SNR}=25,30\mathrm{and}35$ dB, thus yielding three images for each set.

4.1.1. Experimental Setting

For all the comparing methods, the unmixing procedure was performed in a supervised manner, where the ground truth endmembers were utilized. Concerning the initialization of abundance matrix, the result of FCLS was applied in all methods. The parameters in state-of-the-art methods C-SUnSAL, GBM, PPNMM and MLM were basically set according to their original papers. Specifically, a supervised version of MLMp method was considered by simply fixing the endmembers in the algorithm, and the tolerance value was changed to ${10}^{-5}$ from ${10}^{-4}$ as in [18], for the sake of fair comparison. Regarding BNLSU, the maximum iteration number was set to 500, the tolerance value was set to ${10}^{-5}$, and the parameters were set as follows: the penalty parameter in ADMM was set as $\rho =0.05$, and the regularization parameters in Equation (4) were set as ${\lambda }_1=0.001$ for the sparsity constraint and ${\lambda }_2=0.3$ for the smoothness of nonlinearity matrix, following the analysis in [20].

As for the proposed G-MLM algorithm, the maximum iteration number was set to 500, the tolerance value was set to ${10}^{-5}$, and the penalty parameter in ADMM was set as $\rho =0.05$. Four regularization parameters needed to be tuned: $d_{min}^2$, ${\lambda }_1$, ${\lambda }_2$ and ${\lambda }_3$. To keep the comparison fair and also to simplify the process of parameters tuning, we fixed the parameter for sparsity regularization as ${\lambda }_1=0.001$, the same value as in BNLSU.

The threshold parameter $d_{min}^2$ in Equation (6) determines the affinity matrix $\mathit{W}$ computed from the input image, and should be chosen case by case. In essence, this parameter determines the underlying graph structure of the Laplacian regularizers on both abundance matrix and nonlinearity parameter vector, that is, only the pixels connected by nonzero elements in $\mathit{W}$ will mutually influence the parameter estimation. Following the spirit of the work in [36,37], we propose to estimate the parameter $d_{min}^2$ as a function of reconstruction error, by using an empirical formula

$$d_{min}^2=\frac{\theta }{TL}{\parallel \mathit{X}-\mathit{A}{\mathit{S}}_{FCLS}\parallel }_F^2$$

(35)

where the parameter $\theta$ was chosen as $\theta =400$ throughout the experiments, and ${\mathit{S}}_{FCLS}$ is the abundance matrix estimated by FCLS [5]. In this way, the threshold parameter $d_{min}^2$ was chosen as 0.28, 0.17 and 0.13 for DC1 data with $\mathrm{SNR}=25,30,\mathrm{and}35$ dB, respectively, and $d_{min}^2$ was chosen as 0.19, 0.08 and 0.05 for DC2 data with $\mathrm{SNR}=25,30,\mathrm{and}35$ dB, respectively. Of particular note is that the dataset with lower SNR value generally prefers more important threshold value. It is reasonable, since severer noise affects the measure of the real similarity between spectra pairs more.

We verified the appropriateness of the parameter $d_{min}^2$ given by Equation (35), by taking DC1 data at $\mathrm{SNR}=30$ dB as an example. For these data, the affinity matrix $\mathit{W}$ obtained by $d_{min}^2=0.17$ is shown in Figure 2, where a re-arranged matrix is shown for illustration purpose, according to Ammanouil et al. [24]. That is, the pixels of the five squares in the first row are placed first, then the pixels of the five squares in the second row appear, and so on, and the background pixels are placed last. We observe that there are 20 small, 1 medium and 1 big white blocks in $\mathit{W}$. Each white block is a sub-matrix of $\mathit{W}$ whose elements are all close to 1, signifying the associated pixels are mutually similar. This is just consistent with the actual case of data generation. To be precise, the 20 small blocks correspond to the 20 different squares of pixels in the first four rows, each containing $5\times 5=25$ identical pixels; the one medium white block corresponds to an assemblage $25\times 5=125$ identical pixels from five squares in the last row; and the one big white block represents the sub-affinity-matrix for all the remaining identical background pixels. As a result, the graph structure of these data is well-revealed by $\mathit{W}$ obtained by $d_{min}^2=0.17$. In Figure 3, we also compare the affinity matrices $\mathit{W}$ produced by $d_{min}^2=0.20$ and $d_{min}^2=0.28$ given by Equation (35), respectively. It is obvious that, with $d_{min}^2=0.28$, the graph structure is better represented, with the resulting matrix $\mathit{W}$ being more consistent to the actual one. Based on above discussions, the effectiveness of Equation (35) in estimating $d_{min}^2$ is verified to some extent.

The remaining two regularization parameters to be tuned are ${\lambda }_2$ and ${\lambda }_3$, which determine the importances of the graph-based abundance and nonlinearity parameter regularization terms, respectively. To further simplify the parameter tuning process, we applied ${\lambda }_3=\frac{1}{2}{\lambda }_2$ in this experiment, considering that the abundance matrix $\mathit{S}$ has more dimensions than the nonlinear parameter $\mathit{P}$. The parameter ${\lambda }_2$ was then selected by applying the candidate value set $\{0.0001,0.0005,0.001,0.005,0.01,0.05,0.1,0.3,0.5,0.7,0.9,1,1.5,2,2.5,3,3.5,4,4.5,5\}$, on the image with $\mathrm{SNR}=30$ dB for each set. The influence of ${\lambda }_2$ on unmixing performance in terms of RMSE is plotted in Figure 4a,b, where good results are obtained within the interval $[3,4.5]$ and $[0.01,0.1]$ for DC1 and DC2, respectively. Thus, we fixed $({\lambda }_2,{\lambda }_3)=(4,2)$ for DC1, and $({\lambda }_2,{\lambda }_3)=(0.05,0.025)$ for DC2.

4.1.2. Results Analysis

On images in DC1 and DC2, the unmixing performance was evaluated by the RMSE of abundance, as defined in Equation (32). The results obtained by applying all the comparing algorithms on DC1 and DC2 images at different noise levels are compared in

Table 2. The averaged RMSE values over five runs using all the comparing methods, on three DC1 datasets.

Methods	SNR 25 dB	SNR 30 dB	SNR 35 dB
G-MLM	0.0049 ($d_{min}^2=0.28$)	0.0015 ($d_{min}^2=0.17$)	0.0006 ($d_{min}^2=0.13$)
BNLSU	0.0479	0.0206	0.0113
MLMp	0.0194	0.0107	0.0061
MLM	0.0114	0.0065	0.0036
PPNMM	0.0375	0.0365	0.0360
GBM	0.0516	0.0509	0.0506
FCLS	0.0507	0.0499	0.0497

and

Table 3. The averaged RMSE values over five runs by all the comparing methods, on three DC2 datasets.

Methods	SNR 25 dB	SNR 30 dB	SNR 35 dB
G-MLM	0.0302 ($d_{min}^2=0.19$)	0.0285 ($d_{min}^2=0.08$)	0.0282 ($d_{min}^2=0.05$)
BNLSU	0.0504	0.0306	0.0187
MLMp	0.0321	0.0297	0.0290
MLM	0.0311	0.0297	0.0292
PPNMM	0.0128	0.0075	0.0045
GBM	0.0389	0.0381	0.0378
FCLS	0.0388	0.0380	0.0377

, respectively. As observed in Table 2 for DC1, for all the three noise levels, the proposed G-MLM always outperforms its counterparts, by yielding the smallest RMSE value. This result is not surprising, as the proposed G-MLM takes advantage of the graph information hidden in the dataset, which is more consistent with the real situation. Moreover, as the value of SNR increases, the value of RMSE decreases prominently for every method. Namely, high noise level with low SNR deteriorates the performances of all the unmixing strategies. Generally, the MLM-based methods provide better unmixing results than the methods based on linear or bilinear models, on DC1.

In Table 3 for DC2, as the pixels were generated by PPNMM, it is not surprising to observe that the PPNMM unmixing method outperforms all the comparing methods. The proposed G-MLM is able to provide good unmixing results that are only second to those of PPNMM, on DC2 images with $\mathrm{SNR}=25\mathrm{dB}$ and $\mathrm{SNR}=30\mathrm{dB}$. For the image with $\mathrm{SNR}=35\mathrm{dB}$, it is the BNLSU method that yields the second best RMSE.

The estimated abundance maps by using different methods are shown in Figure 5 for DC1 image with SNR = 30, and in Figure 6 for DC2 image with SNR = 30. It is not difficult to find that the abundance maps obtained by the proposed G-MLM are the closest to the ground truth, when compared to all the other methods. For DC1 image with SNR = 30, the corresponding maps of P obtained by the MLM-based unmixing methods are shown in Figure 7, where the proposed G-MLM leads to the P map most consistent to the ground turth.

To study the convergence of the proposed ADMM algorithm for G-MLM, we traced the value of the objective function in Equation (9), the reconstruction error, the primal residual and the dual residual along with the number of iterations, as given in Figure 8. It is observed that all these values drop rapidly for the first 50 iterations, and then nearly decline to zeros at a slower rate. It demonstrates the good convergence of the proposed ADMM algorithm, according to [6,29].

Lastly, the executing time by different methods are compared in

Table 4. Computational time measured by different algorithms on DC1 with SNR = 30 dB.

Methods	Time [s]
G-MLM	445.18
BNLSU	574.31
MLMp	69.43
MLM	118.88
PPNMM	235.50
GBM	34.25
FCLS	1.93

. The original G-MLM is shown to be time-consuming. To alleviate this issue, we examine ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ and ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ algorithms. The comparison of these three methods are reported in

Table 5. Comparison of RMSE and computational time of G-MLM, ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ and ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ on DC1 with SNR = 30 dB. The preprocessing time for ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ and ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ is omitted in the table.

Methods	RMSE	Time [s]
G-MLM	0.0015	445.18
G-MLMsub	0.0015	402.29
G-MLMsuper	0.0027	152.42

, where the number of clusters in ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ is set as $k=10$ and the number of superpixels in ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ is set as $m=120$. As observed, the use of ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ algorithm slightly reduces the computational time of G-MLM, almost without deteriorating the unmixing performance. The ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ algorithm can greatly save the computational time, while providing acceptable RMSE that is still smaller than that of the comparing methods in Table 2. It is noteworthy that the overriding merit of ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ is the practicality for addressing large-scale images, especially in the cases where both G-MLM and ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ may not be applicable.

4.2. Experiments on Urban Dataset

The Urban dataset is widely-investigated in hyperspectral unmixing research [9,38]. The original image has $307\times 307$ pixels with each pixel corresponding to an area of $2\times 2$ m${}^2$ area. It is measured by 210 wavelengths ranging from 400 nm to 2500 nm. After removing the noisy bands (1–4, 76, 87, 101–111, 136–153, and 198–210) that are contaminated by water vapor density and atmosphere, the remaining 162 bands of high SNR were considered. This field is known to be mainly composed by four materials: roof, tree, asphalt and grass. In our experiment, the four endmembers were first extracted by the VCA algorithm, as shown in Figure 9.

Because the Urban image is of large scale, neither G-MLM nor ${\mathrm{G}\text{-}\mathrm{MLM}}_{sub}$ can be directly applied on our computer. Therefore, the ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ scheme was adopted. We speciied the selection of parameters in ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$ for Urban image. The original image was firstly divided into $m=120$ superpixels, as shown in Figure 10. The value of $d_{min}^2$ in Equation (6) was estimated to be 0.02, according to Equation (35). As done for synthetic images, the sparsity regularization parameter was set as ${\lambda }_1=0.001$, and the relationship with ${\lambda }_3=\frac{1}{2}{\lambda }_2$ was maintained. For Urban image, the value of ${\lambda }_2$ was selected from the set $\{0.0001,0.0005,0.001,0.005,0.01,0.05,0.1,0.5,1,2,3,4,5\}$. For illustration purposes, we plot the changes of RE along with ${\lambda }_2$ in Figure 11. It is observed that the RE value is relatively small and stable for small values of ${\lambda }_2$, and starts to increase when ${\lambda }_2$ exceeds 0.001. Therefore, we chose $({\lambda }_2,{\lambda }_3)=(0.001,0.0005)$ for Urban image.

The executing time using all the comparing methods is given in

Table 6. Comparison of RE and computational time of ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$, BNLSU, MLMp, MLM, PPNMM, GBM, and FCLS on Urabn image.

Methods	RE	Time [s]
${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$	0.0053	1904.89
BNLSU	0.0019	5443.02
MLMp	0.0053	26.77
MLM	0.0053	1075.40
PPNMM	0.0052	591.42
GBM	0.0056	72.37
FCLS	0.0056	19.09

. A comparison of RE is also reported in the same table, where BNLSU provides a value much smaller than other unmixing approaches. It is reasonable, as this method defines the nonlinearity parameter in each band for a given pixel, thus yielding a model with many more parameters to fit the reconstruction error, compared to other methods. However, it should be noticed that RE is only adopted as a reference of the unmixing performance. It measures the averaged fitting level of a mixture model, e.g., LMM, GBM or MLM, to the observed spectra, but cannot truly reflect the unmixing quality in terms of abundance, especially for real hyperspectral images where the mechanism of mixture is complicated and unclear [10]. The abundance maps obtained by all comparing algorithms are shown in Figure 12. The abundance maps by different algorithms visually show little differences. The estimated maps of P obtained with ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$, MLM and MLMp are shown in Figure 13. All three maps correspond well with the areas in the image where nonlinear effects are expected to be significant, e.g., the intersection area between roof and tree, and between grass and asphalt. As the graph structure of data is considered in the proposed ${\mathrm{G}\text{-}\mathrm{MLM}}_{super}$, the corresponding abundance maps and P map exhibit more spatial consistency and and smoothness compared with other algorithms.

5. Conclusions

This paper proposes a graph regularized hyperspectral unmixing method based on the multilinear mixing model. By taking advantage of the graph structure embedded in the data, Laplacian graph ragularizers are introduced to regularize both the abundance matrix and the nonlinear probability parameters. The sparse constraint is also enforced on the abundance matrix. The resulting optimization problem is solved by using ADMM algorithm, and the superpixel construction scheme is applied to reduce the complexity. Comparative studies on two series of synthetic data and a real hyperspectral image verified the advantage of the proposed algorithm in terms of unmixing results. Future works include jointly exploiting the graph and spatial structures to further augment the G-MLM model, in order to accurately characterize abundances and nonlinear probability parameters. Extending the nonlinear parameters in G-MLM model to a band-wise form is also worth our investigation.