# Image Processing and Machine Learning for Hyperspectral Unmixing: An Overview and the HySUPP Python Package

Behnood Rasti, Alexandre Zouaoui, Julien Mairal, Jocelyn Chanussot

► To cite this version:

Behnood Rasti, Alexandre Zouaoui, Julien Mairal, Jocelyn Chanussot. Image Processing and Machine Learning for Hyperspectral Unmixing: An Overview and the HySUPP Python Package. IEEE Transactions on Geoscience and Remote Sensing, 2024, pp.1-31. 10.1109/TGRS.2024.3393570. hal-04180307v3

HAL Id: hal-04180307

https://hal.science/hal-04180307v3

Submitted on 24 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Image Processing and Machine Learning for Hyperspectral Unmixing: An Overview and the HySUPP Python Package

Behnood Rasti, Alexandre Zouaoui, Julien Mairal, Jocelyn Chanussot

## ► To cite this version:

Behnood Rasti, Alexandre Zouaoui, Julien Mairal, Jocelyn Chanussot. Image Processing and Machine Learning for Hyperspectral Unmixing: An Overview and the HySUPP Python Package. IEEE Transactions on Geoscience and Remote Sensing, 2024, 62 (5517631), pp.1-31. 10.1109/TGRS.2024.3393570. [hal-04180307v3](https://hal.science/hal-04180307v3)

HAL Id: hal-04180307

https://hal.science/hal-04180307v3

Submitted on 24 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Image Processing and Machine Learning for Hyperspectral Unmixing: An Overview and the HySUPP Python Package

Behnood Rasti, *Senior Member, IEEE*, Alexandre Zouaoui, *Student Member, IEEE*, Julien Mairal, *Senior Member, IEEE* and Jocelyn Chanussot, *Fellow, IEEE*

***Abstract***—Spectral pixels are often a mixture of the pure spectra of the materials, called endmembers, due to the low spatial resolution of hyperspectral sensors, double scattering, and intimate mixtures of materials in the scenes. Unmixing estimates the fractional abundances of the endmembers within the pixel. Depending on the prior knowledge of endmembers, linear unmixing can be divided into three main groups: supervised, semi-supervised, and unsupervised (blind) linear unmixing. Advances in image processing and machine learning substantially affected unmixing. This paper provides an overview of advanced and conventional unmixing approaches. Additionally, we draw a critical comparison between advanced and conventional techniques from the three categories. We compare the performance of the unmixing techniques on three simulated and one real dataset. The experimental results reveal the advantages of different unmixing categories for different unmixing scenarios. Moreover, we provide an open-source Python-based package available at https://github.com/BehnoodRasti/HySUPP to reproduce the results.

***Index Terms***—Hyperspectral, unmixing, endmember extraction, abundance estimation, linear mixture, machine learning, deep learning, optimization.

## I. INTRODUCTION

SPECTRAL unmixing is a crucial processing technique in hyperspectral remote sensing that can play a vital role in various fields such as mineral exploration, agriculture and crop monitoring, environmental monitoring, urban planning, remote sensing of planetary surfaces, pollution monitoring, medical imaging, water quality assessment, etc. The ability to separate and identify different materials in an image is made possible by the contiguous spectra captured by hyperspectral sensors. Through the use of endmembers, which are the unique spectral signatures of materials, unmixing algorithms can decompose the mixed spectral data into its constituent parts. However, due to low spatial resolution, multiple scattering, and intimate mixing, the measured spectrum within a pixel is generally a complex mixture of the pure spectra of the constituent materials, making unmixing a challenging task. Fig. 1 demonstrates how the reflectance of a mixed pixel

captured by an optical hyperspectral camera is composed of two endmembers within that pixel. In hyperspectral remote sensing, a mixing model represents the observed spectral pixel as a function of the endmembers and their corresponding fractional abundances within the pixel's area. Unmixing is the process of estimating the fractional abundances, either by estimating or extracting the endmembers or by relying on a library of endmembers. It may also involve determining the number of endmembers present. The mixing model is either linear or nonlinear, depending on the interaction of the incident light and the materials in the scene or sample.

In linear unmixing, the endmembers are assumed to be linearly mixed, which is valid when each light ray interacts with only one material before reaching the sensor, as shown in Figure 2 (a). This assumption is common in Earth observation applications where macroscopic problems exist. The spatial resolution of the sensor plays a crucial role in macroscopic scenarios, where a pixel may contain multiple materials. In such cases, the pixel is a mixture of different materials, such as tree and soil, as shown in Figure 2 (a).

Another assumption is bilinear mixing, which assumes double scattering or that the light ray interacts with two materials before reaching the sensor, as shown in Figure 2 (b). However, in microscopic scenarios where the pure materials are intimately mixed within a pixel, and the light undergoes multiple scattering and reflections by several materials, the linear approximation fails, and nonlinear models that are more complex than bilinear should be used [1], [2].

In macroscopic scenarios, nonlinearity is also induced by noise, atmospheric effects, temporal effects, illumination variations, and materials' intrinsic variation, which may cause spectral variability shown in Fig. 3 [3]. Illumination variations are mainly due to two effects: varying terrain topography, which affects the angles of the incident radiation, and the occlusion of the light source by other objects (leading to shaded areas). The intrinsic variation of spectral signatures of the materials is mainly due to physicochemical changes. For instance, soil's signature might change dramatically by variations in its composition and moisture content. Another example is the leaf's signature changes throughout the temporal seasons. Scale differences can cause spectral variability as well. For instance, a tree contains leaves, branches, fruits, bark, bark, fruits or flowers, etc., leading to the intrinsic variability of the tree spectra at different scales. The atmospheric effect should be corrected, and noise can be removed [4]–[6].

Behnood Rasti (corresponding author) is with the Faculty of Electrical Engineering and Computer Science, Technische Universität Berlin, 10623 Berlin, Germany; behnood.rasti@gmail.com

Alexandre Zouaoui, Jocelyn Chanussot, and Julien Mairal are with Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France.

Alexandre Zouaoui is currently with Data Science Experts, Grenoble, France.

Manuscript received April 19, 2023; revised August 16, 2023.
# Image Processing and Machine Learning for Hyperspectral Unmixing: An Overview and the HySUPP Python Package

Behnood Rasti, *Senior Member, IEEE*, Alexandre Zouaoui, *Student Member, IEEE*, Julien Mairal, *Senior Member, IEEE* and Jocelyn Chanussot, *Fellow, IEEE*

***Abstract***—Spectral pixels are often a mixture of the pure spectra of the materials, called endmembers, due to the low spatial resolution of hyperspectral sensors, double scattering, and intimate mixtures of materials in the scenes. Unmixing estimates the fractional abundances of the endmembers within the pixel. Depending on the prior knowledge of endmembers, linear unmixing can be divided into three main groups: supervised, semi-supervised, and unsupervised (blind) linear unmixing. Advances in image processing and machine learning substantially affected unmixing. This paper provides an overview of advanced and conventional unmixing approaches. Additionally, we draw a critical comparison between advanced and conventional techniques from the three categories. We compare the performance of the unmixing techniques on three simulated and one real dataset. The experimental results reveal the advantages of different unmixing categories for different unmixing scenarios. Moreover, we provide an open-source Python-based package available at https://github.com/BehnoodRasti/HySUPP to reproduce the results.

***Index Terms***—Hyperspectral, unmixing, endmember extraction, abundance estimation, linear mixture, machine learning, deep learning, optimization.

## I. INTRODUCTION

SPECTRAL unmixing is a crucial processing technique in hyperspectral remote sensing that can play a vital role in various fields such as mineral exploration, agriculture and crop monitoring, environmental monitoring, urban planning, remote sensing of planetary surfaces, pollution monitoring, medical imaging, water quality assessment, etc. The ability to separate and identify different materials in an image is made possible by the contiguous spectra captured by hyperspectral sensors. Through the use of endmembers, which are the unique spectral signatures of materials, unmixing algorithms can decompose the mixed spectral data into its constituent parts. However, due to low spatial resolution, multiple scattering, and intimate mixing, the measured spectrum within a pixel is generally a complex mixture of the pure spectra of the constituent materials, making unmixing a challenging task. Fig. 1 demonstrates how the reflectance of a mixed pixel

captured by an optical hyperspectral camera is composed of two endmembers within that pixel. In hyperspectral remote sensing, a mixing model represents the observed spectral pixel as a function of the endmembers and their corresponding fractional abundances within the pixel's area. Unmixing is the process of estimating the fractional abundances, either by estimating or extracting the endmembers or by relying on a library of endmembers. It may also involve determining the number of endmembers present. The mixing model is either linear or nonlinear, depending on the interaction of the incident light and the materials in the scene or sample.

In linear unmixing, the endmembers are assumed to be linearly mixed, which is valid when each light ray interacts with only one material before reaching the sensor, as shown in Figure 2 (a). This assumption is common in Earth observation applications where macroscopic problems exist. The spatial resolution of the sensor plays a crucial role in macroscopic scenarios, where a pixel may contain multiple materials. In such cases, the pixel is a mixture of different materials, such as tree and soil, as shown in Figure 2 (a).

Another assumption is bilinear mixing, which assumes double scattering or that the light ray interacts with two materials before reaching the sensor, as shown in Figure 2 (b). However, in microscopic scenarios where the pure materials are intimately mixed within a pixel, and the light undergoes multiple scattering and reflections by several materials, the linear approximation fails, and nonlinear models that are more complex than bilinear should be used [1], [2].

In macroscopic scenarios, nonlinearity is also induced by noise, atmospheric effects, temporal effects, illumination variations, and materials' intrinsic variation, which may cause spectral variability shown in Fig. 3 [3]. Illumination variations are mainly due to two effects: varying terrain topography, which affects the angles of the incident radiation, and the occlusion of the light source by other objects (leading to shaded areas). The intrinsic variation of spectral signatures of the materials is mainly due to physicochemical changes. For instance, soil's signature might change dramatically by variations in its composition and moisture content. Another example is the leaf's signature changes throughout the temporal seasons. Scale differences can cause spectral variability as well. For instance, a tree contains leaves, branches, fruits, bark, bark, fruits or flowers, etc., leading to the intrinsic variability of the tree spectra at different scales. The atmospheric effect should be corrected, and noise can be removed [4]–[6].

Behnood Rasti (corresponding author) is with the Faculty of Electrical Engineering and Computer Science, Technische Universität Berlin, 10623 Berlin, Germany; behnood.rasti@gmail.com

Alexandre Zouaoui, Jocelyn Chanussot, and Julien Mairal are with Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France.

Alexandre Zouaoui is currently with Data Science Experts, Grenoble, France.

Manuscript received April 19, 2023; revised August 16, 2023.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparisons of mixed and pure pixels in hyperspectral data.</figcaption></figure>

The unmixing process is subjective and depends on how the endmembers are defined. For example, when detecting buildings, the reflectance of a brick could be considered an endmember. However, bricks are composed of several materials, such as clay, sand, and concrete, making it challenging to define the endmembers of interest for detecting specific materials or ratios of different materials inside bricks.

In the former, linear unmixing may be sufficient to solve the problem. However, for intimate mixtures and the latter, nonlinear models must be used. For example, the Hapke model [7], [8] suggests that mixing occurs at the albedo level rather than the reflectance level. Therefore, more complex nonlinear models are necessary for such cases to accurately unmix the spectral data. It is worth mentioning that, unlike the common assumption, the proportion of an endmember within a pixel or in a scene is not the percentage of the endmember material within that pixel or in that scene [1]. According to Hapke [7], [8], the proportion of an endmember in a linear mixture indicates the relative area of that endmember. Therefore, reflectance is not a linear mixture of the mass or cross-sectional area of the endmember materials. Despite this, linear unmixing has shown significant value in remote sensing applications in the past decades [1]. However, when associating an unmixing problem with an application, we should be aware of the simplifications of the model we use.

Assuming three materials in a scene correspond to three endmembers ($e_1$, $e_2$, and $e_3$), the observed data can be projected into a 2D subspace where all the data points are enclosed in the data simplex formed by the convex hull of endmembers (see Fig. 11). In contrast, in a nonlinear scenario, the simplex does not contain all the data points. This paper mainly considers linear EO applications and macroscopic problems and focuses on linear mixing models. However, we briefly discuss the bilinear and nonlinear unmixing approaches.

### A. Key trends in the hyperspectral unmixing realm

Hyperspectral unmixing has undergone various developmental stages, as illustrated in figure 4. In its early years, from 1992 to 2008, this research topic received minimal attention, with an annual publication count of no more than 18. However, starting from 2009 and continuing until 2016, there was a steady increase in interest, with the number of publications even surpassing 200 in 2016. Since 2017, hyperspectral unmixing has remained a prominent and dynamic area of research, driven by the emergence of deep learning techniques for unmixing. In 2022, this field reached its fourth most prolific year on record, with 164 publications.

During the period when hyperspectral unmixing gained popularity, the primary scientific programming language of choice was MATLAB, as depicted in figure 5. However, in recent times, there has been a steady increase in interest and adoption of Python, while MATLAB's popularity has waned. This shift can be attributed, in part, to the growing trend of open science, where open-sourcing scientific code has become widely embraced. Consequently, the present scenario suggests that developing and releasing an open-source unmixing package in Python is more advantageous than doing so in MATLAB, even though it may involve the translation of existing unmixing code from MATLAB to Python.

### B. Existing Overview Publications and Unmixing Packages

An early survey on hyperspectral unmixing was given in [9], which discusses basic geometrical and statistical methods. In [1], linear model-based unmixing techniques were divided into three categories: geometrical, statistical, and sparse regression-based approaches. A Matlab toolbox is available at https://openremotesensing.net/. However, the toolbox is incomplete, and some methods, such as dependent component analysis (DECA), are missing. In recent years, deep learning and neural networks have become state-of-the-art for many tasks in machine learning and image processing. Consequently, many unmixing approaches were proposed based on shallow and deep neural networks. A comparison of autoencoder-based networks was drawn in [10]. The authors discussed autoencoder-based architectures divided into five categories, i.e., Sparse Nonnegative Autoencoders, Variational Autoencoders, Adversarial Autoencoders, Denoising Autoencoders, and Convolutional Autoencoders. They further discuss the choice of different modules, such as different activation functions or loss functions, and they compared shallow networks to deep ones. They also provide a TensorFlow-based Python package that is available on GitHub. However, the package is limited to blind unmixing approaches based on autoencoders. It does not discuss or compare the supervised, semisupervised, and more conventional blind unmixing methods.

In [11], some model-based and neural network-based unmixing approaches were explained but without experimental comparisons. A list of resources for the approaches was given.

A survey on endmember variability in Spectral Mixture Analysis (SMA) was given in [12]. In [13], an overview of unmixing methods that address endmember variability was given. A comprehensive overview of the unmixing methods
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparisons of mixed and pure pixels in hyperspectral data.</figcaption></figure>

The unmixing process is subjective and depends on how the endmembers are defined. For example, when detecting buildings, the reflectance of a brick could be considered an endmember. However, bricks are composed of several materials, such as clay, sand, and concrete, making it challenging to define the endmembers of interest for detecting specific materials or ratios of different materials inside bricks.

In the former, linear unmixing may be sufficient to solve the problem. However, for intimate mixtures and the latter, nonlinear models must be used. For example, the Hapke model [7], [8] suggests that mixing occurs at the albedo level rather than the reflectance level. Therefore, more complex nonlinear models are necessary for such cases to accurately unmix the spectral data. It is worth mentioning that, unlike the common assumption, the proportion of an endmember within a pixel or in a scene is not the percentage of the endmember material within that pixel or in that scene [1]. According to Hapke [7], [8], the proportion of an endmember in a linear mixture indicates the relative area of that endmember. Therefore, reflectance is not a linear mixture of the mass or cross-sectional area of the endmember materials. Despite this, linear unmixing has shown significant value in remote sensing applications in the past decades [1]. However, when associating an unmixing problem with an application, we should be aware of the simplifications of the model we use.

Assuming three materials in a scene correspond to three endmembers ($e_1$, $e_2$, and $e_3$), the observed data can be projected into a 2D subspace where all the data points are enclosed in the data simplex formed by the convex hull of endmembers (see Fig. 11). In contrast, in a nonlinear scenario, the simplex does not contain all the data points. This paper mainly considers linear EO applications and macroscopic problems and focuses on linear mixing models. However, we briefly discuss the bilinear and nonlinear unmixing approaches.

### A. Key trends in the hyperspectral unmixing realm

Hyperspectral unmixing has undergone various developmental stages, as illustrated in figure 4. In its early years, from 1992 to 2008, this research topic received minimal attention, with an annual publication count of no more than 18. However, starting from 2009 and continuing until 2016, there was a steady increase in interest, with the number of publications even surpassing 200 in 2016. Since 2017, hyperspectral unmixing has remained a prominent and dynamic area of research, driven by the emergence of deep learning techniques for unmixing. In 2022, this field reached its fourth most prolific year on record, with 164 publications.

During the period when hyperspectral unmixing gained popularity, the primary scientific programming language of choice was MATLAB, as depicted in figure 5. However, in recent times, there has been a steady increase in interest and adoption of Python, while MATLAB's popularity has waned. This shift can be attributed, in part, to the growing trend of open science, where open-sourcing scientific code has become widely embraced. Consequently, the present scenario suggests that developing and releasing an open-source unmixing package in Python is more advantageous than doing so in MATLAB, even though it may involve the translation of existing unmixing code from MATLAB to Python.

### B. Existing Overview Publications and Unmixing Packages

An early survey on hyperspectral unmixing was given in [9], which discusses basic geometrical and statistical methods. In [1], linear model-based unmixing techniques were divided into three categories: geometrical, statistical, and sparse regression-based approaches. A Matlab toolbox is available at https://openremotesensing.net/. However, the toolbox is incomplete, and some methods, such as dependent component analysis (DECA), are missing. In recent years, deep learning and neural networks have become state-of-the-art for many tasks in machine learning and image processing. Consequently, many unmixing approaches were proposed based on shallow and deep neural networks. A comparison of autoencoder-based networks was drawn in [10]. The authors discussed autoencoder-based architectures divided into five categories, i.e., Sparse Nonnegative Autoencoders, Variational Autoencoders, Adversarial Autoencoders, Denoising Autoencoders, and Convolutional Autoencoders. They further discuss the choice of different modules, such as different activation functions or loss functions, and they compared shallow networks to deep ones. They also provide a TensorFlow-based Python package that is available on GitHub. However, the package is limited to blind unmixing approaches based on autoencoders. It does not discuss or compare the supervised, semisupervised, and more conventional blind unmixing methods.

In [11], some model-based and neural network-based unmixing approaches were explained but without experimental comparisons. A list of resources for the approaches was given.

A survey on endmember variability in Spectral Mixture Analysis (SMA) was given in [12]. In [13], an overview of unmixing methods that address endmember variability was given. A comprehensive overview of the unmixing methods
<figure><img src="image_3.png" /><figcaption>Fig. 2. Macroscopic versus microscopic assumption. This figure illustrates three major assumptions in hyperspectral imagery. (a) The linear assumption is when the light interacts with the materials only once. (b) The bilinear assumption when the light interacts with a maximum of two materials (c) intimate mixture when the light interacts with more than two materials.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Noise, atmospheric effects, illumination variations (caused by terrain topography, occlusion of the light), and intrinsic variations of materials (e.g., soil signature might change dramatically by variations in its composition and moisture content) cause spectral variability.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Publications over time based on IEEE Xplore keyword search tool using “Hyperspectral Unmixing” as input.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Worldwide interest in scientific programming languages over time according to Google Trends in the “Science” category.</figcaption></figure>

that address spectral variability was recently provided in [3], and a list of Matlab codes was also given. In [14], an overview of endmember extraction approaches was given. Review papers on hyperspectral remote sensing data analysis briefly discussed the unmixing methods [15], [16]. In [17], an overview of tensor-based unmixing models was given. In [18], a survey of nonlinear unmixing methods is given.

There are other open-source tools such as HyperMix [19], Spectral Python (SPy) (https://www.spectralpython.net/), Spectral Library Tool (https://spectral-libraries.readthedocs.io/en/latest/), PySptools (https://pypi.org/project/pysptools/), that include basic algorithms for estimation of the number of endmembers, endmember extraction, abundance estimation, and some library tools, and library-based methods. Thus, there is a need for a comprehensive package that covers the methodologies across different unmixing categories and contains state-of-the-art image processing and machine learning techniques.

### C. Significance of Contribution

This paper aims at the followings:

* Providing researchers with a comprehensive yet technical overview of all the essential topics regarding linear unmixing techniques.

* Categorizing the unmixing approaches considering prior knowledge available about endmembers. Linear unmixing can be divided into three main categories: supervised, semi-supervised (library-based), and unsupervised (blind) unmixing.

* Compare the unmixing methods in terms of prior knowledge of the endmembers and draw conclusions that can help researchers to select an appropriate unmixing method to tackle real-world challenges. We compare conventional and deep learning-based unmixing approaches in those categories for three simulated and two real-world datasets. For the simulated datasets, we consider three scenarios: a simple, pure pixel dataset, a dataset with spectral variations, and a challenging dataset with no pure pixel. Such comparisons provide insight to researchers into which category to use for their application. Additionally, the comparisons reveal the drawbacks of the categories, which motivate the developers to investigate new ideas to address them.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Macroscopic versus microscopic assumption. This figure illustrates three major assumptions in hyperspectral imagery. (a) The linear assumption is when the light interacts with the materials only once. (b) The bilinear assumption when the light interacts with a maximum of two materials (c) intimate mixture when the light interacts with more than two materials.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Noise, atmospheric effects, illumination variations (caused by terrain topography, occlusion of the light), and intrinsic variations of materials (e.g., soil signature might change dramatically by variations in its composition and moisture content) cause spectral variability.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Publications over time based on IEEE Xplore keyword search tool using “Hyperspectral Unmixing” as input.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Worldwide interest in scientific programming languages over time according to Google Trends in the “Science” category.</figcaption></figure>

that address spectral variability was recently provided in [3], and a list of Matlab codes was also given. In [14], an overview of endmember extraction approaches was given. Review papers on hyperspectral remote sensing data analysis briefly discussed the unmixing methods [15], [16]. In [17], an overview of tensor-based unmixing models was given. In [18], a survey of nonlinear unmixing methods is given.

There are other open-source tools such as HyperMix [19], Spectral Python (SPy) (https://www.spectralpython.net/), Spectral Library Tool (https://spectral-libraries.readthedocs.io/en/latest/), PySptools (https://pypi.org/project/pysptools/), that include basic algorithms for estimation of the number of endmembers, endmember extraction, abundance estimation, and some library tools, and library-based methods. Thus, there is a need for a comprehensive package that covers the methodologies across different unmixing categories and contains state-of-the-art image processing and machine learning techniques.

### C. Significance of Contribution

This paper aims at the followings:

* Providing researchers with a comprehensive yet technical overview of all the essential topics regarding linear unmixing techniques.

* Categorizing the unmixing approaches considering prior knowledge available about endmembers. Linear unmixing can be divided into three main categories: supervised, semi-supervised (library-based), and unsupervised (blind) unmixing.

* Compare the unmixing methods in terms of prior knowledge of the endmembers and draw conclusions that can help researchers to select an appropriate unmixing method to tackle real-world challenges. We compare conventional and deep learning-based unmixing approaches in those categories for three simulated and two real-world datasets. For the simulated datasets, we consider three scenarios: a simple, pure pixel dataset, a dataset with spectral variations, and a challenging dataset with no pure pixel. Such comparisons provide insight to researchers into which category to use for their application. Additionally, the comparisons reveal the drawbacks of the categories, which motivate the developers to investigate new ideas to address them.
* We provide an open-source HyperSpectral Unmixing Python Package (HySUPP). HySUPP is the first open-source python-based hyperspectral unmixing package to include supervised, semi-supervised, and blind unmixing methods. The package will benefit the geoscience and remote sensing community, including researchers, developers, lecturers, and students. The package installation is straightforward since HySUPP relies on a few dependencies. In addition, all the methods can be run using a single command line instruction.

## II. HYSUPP: HYPERSPECTRAL UNMIXING PYTHON PACKAGE

HySUPP exhibits a list of highly desirable properties summarized as follows: i) completeness, ii) reproducibility, iii) extensibility and iv) homogeneity. It implements common best practices and enables simple benchmarking of unmixing techniques thanks to user-friendly command line instructions.

### A. Features

a) **Completeness:** As a practitioner, the ability to explore and experiment with different unmixing techniques is crucial since no single approach can consistently outperform others in all unmixing scenarios. Thus, ensuring the completeness of our toolbox becomes essential. Our toolbox is designed to cover all three types of unmixing - supervised, semi-supervised, and unsupervised - while striving to be as representative as possible of the various unmixing approaches, although aiming for exhaustive inclusion would be impractical. Currently, HySUPP provides access to a diverse set of 20 different unmixing methods, including 6 supervised, 6 semi-supervised, and 8 unsupervised techniques.

b) **Reproducibility:** Ensuring experimental reproducibility is crucial when exploring various unmixing techniques, as it guarantees the robustness of the conclusions drawn by users. In line with this principle, our toolbox offers the ability to seed experiments, facilitating reproducibility through repeatable noise generation. Additionally, HySUPP automatically saves estimates outputs, providing users with a convenient way to review and compare results, thus enhancing the reliability of their research findings.

c) **Extensibility:** HySUPP's architecture is designed to support the effortless integration of new methods, providing a platform for future advancements in hyperspectral unmixing. Leveraging configuration files, we empower researchers to conduct experiments with ease and flexibility, enabling them to effortlessly incorporate their own models into the toolbox.

d) **Homogeneity:** The uniformity of HySUPP's codebase ensures consistency across different components, enhancing its overall usability and reliability. More specifically, we use common methods for models and common attributes for datasets.

e) **Best practices:** We offer a straightforward yet potent method to monitor the objective function of each approach using Python's `tqdm` library. Furthermore, our pipeline incorporates precise endmember auto-alignment through the utilization of the `munkres` algorithm, provided ground truth

abundance maps are accessible. This enhancement ensures the accurate computation of unmixing performance. By establishing these best practices, we streamline the implementation process and foster a collaborative environment where researchers can easily build upon existing work and share their contributions effectively.

f) **Benchmarking:** Owing to its rich set of features, HySUPP enables users to easily benchmark all implemented methods on their dataset of choice. Our toolbox currently provides 3 synthetic datasets corresponding to different unmixing scenarios. Moreover, we incorporate 4 distinct metrics to thoroughly evaluate unmixing accuracy across methods. Finally, leveraging `mlxp` [20] results query tool, we empower users to analyze and visualize their results in an appealing and informative manner.

### B. Example

The following command line instruction provides an example on how to run a semi-supervised unmixing technique, `SUnCNN` [21], on the `DC1` dataset using an optional custom value for the signal-to-noise ratio (SNR):

```
$ python unmixing.py mode=semi data=DC1 model=SUnCNN noise.SNR=30
```

Table I lists unmixing methods included in HySUPP with their corresponding dependencies. We highlighted our main contributions and the link to the original implementations.

## III. LINEAR UNMIXING

For frequent Earth observation, hyperspectral sensors are installed on satellites. The light passes through the atmosphere before reaching the target and then passes through the atmosphere to reach the sensor. The atmosphere absorbs and scatters light, and therefore, atmospheric corrections should be applied in addition to converting radiance to reflectance. Arguably, atmospheric corrections are unnecessary for airborne, drone-borne, and close-range imaging.

Fig. 6 (a) simplifies the concept of sensing using a satellite hyperspectral sensor. As can be seen, the sensor captures a pixel that contains three materials, i.e., Water, Tree, and Soil. At the sensor, irradiance is corrected to the reflectance in Fig. 6 (b). Note that this reflectance is between zero and one. Fig. 6 (c) shows the concept of linear unmixing, which demonstrates that the measured pixel contains 20% water (blue endmember), 50% soil (yellow endmember), and 30% tree (green endmember). Arguably, these portions cannot be negative, known as abundance nonnegativity constraint (ANC), and the summation of those portions is 100%, known as abundance sum to one constraint (ASC). Additionally, the endmembers are reflectance and cannot be negative (indeed, the endmembers have to be bounded between 0 and 1), which led to the endmember nonnegativity constraint. Here, we discuss two widely used linear mixing models, i.e., the low-rank linear mixing model (so-called LMM) and the sparse and redundant linear model. The former was mainly used for supervised and blind unmixing and the latter for sparse unmixing. Some other linear models have been used in the
* We provide an open-source HyperSpectral Unmixing Python Package (HySUPP). HySUPP is the first open-source python-based hyperspectral unmixing package to include supervised, semi-supervised, and blind unmixing methods. The package will benefit the geoscience and remote sensing community, including researchers, developers, lecturers, and students. The package installation is straightforward since HySUPP relies on a few dependencies. In addition, all the methods can be run using a single command line instruction.

## II. HYSUPP: HYPERSPECTRAL UNMIXING PYTHON PACKAGE

HySUPP exhibits a list of highly desirable properties summarized as follows: i) completeness, ii) reproducibility, iii) extensibility and iv) homogeneity. It implements common best practices and enables simple benchmarking of unmixing techniques thanks to user-friendly command line instructions.

### A. Features

a) **Completeness:** As a practitioner, the ability to explore and experiment with different unmixing techniques is crucial since no single approach can consistently outperform others in all unmixing scenarios. Thus, ensuring the completeness of our toolbox becomes essential. Our toolbox is designed to cover all three types of unmixing - supervised, semi-supervised, and unsupervised - while striving to be as representative as possible of the various unmixing approaches, although aiming for exhaustive inclusion would be impractical. Currently, HySUPP provides access to a diverse set of 20 different unmixing methods, including 6 supervised, 6 semi-supervised, and 8 unsupervised techniques.

b) **Reproducibility:** Ensuring experimental reproducibility is crucial when exploring various unmixing techniques, as it guarantees the robustness of the conclusions drawn by users. In line with this principle, our toolbox offers the ability to seed experiments, facilitating reproducibility through repeatable noise generation. Additionally, HySUPP automatically saves estimates outputs, providing users with a convenient way to review and compare results, thus enhancing the reliability of their research findings.

c) **Extensibility:** HySUPP's architecture is designed to support the effortless integration of new methods, providing a platform for future advancements in hyperspectral unmixing. Leveraging configuration files, we empower researchers to conduct experiments with ease and flexibility, enabling them to effortlessly incorporate their own models into the toolbox.

d) **Homogeneity:** The uniformity of HySUPP's codebase ensures consistency across different components, enhancing its overall usability and reliability. More specifically, we use common methods for models and common attributes for datasets.

e) **Best practices:** We offer a straightforward yet potent method to monitor the objective function of each approach using Python's `tqdm` library. Furthermore, our pipeline incorporates precise endmember auto-alignment through the utilization of the `munkres` algorithm, provided ground truth

abundance maps are accessible. This enhancement ensures the accurate computation of unmixing performance. By establishing these best practices, we streamline the implementation process and foster a collaborative environment where researchers can easily build upon existing work and share their contributions effectively.

f) **Benchmarking:** Owing to its rich set of features, HySUPP enables users to easily benchmark all implemented methods on their dataset of choice. Our toolbox currently provides 3 synthetic datasets corresponding to different unmixing scenarios. Moreover, we incorporate 4 distinct metrics to thoroughly evaluate unmixing accuracy across methods. Finally, leveraging `mlxp` [20] results query tool, we empower users to analyze and visualize their results in an appealing and informative manner.

### B. Example

The following command line instruction provides an example on how to run a semi-supervised unmixing technique, `SUnCNN` [21], on the `DC1` dataset using an optional custom value for the signal-to-noise ratio (SNR):

```
$ python unmixing.py mode=semi data=DC1 model=SUnCNN noise.SNR=30
```

Table I lists unmixing methods included in HySUPP with their corresponding dependencies. We highlighted our main contributions and the link to the original implementations.

## III. LINEAR UNMIXING

For frequent Earth observation, hyperspectral sensors are installed on satellites. The light passes through the atmosphere before reaching the target and then passes through the atmosphere to reach the sensor. The atmosphere absorbs and scatters light, and therefore, atmospheric corrections should be applied in addition to converting radiance to reflectance. Arguably, atmospheric corrections are unnecessary for airborne, drone-borne, and close-range imaging.

Fig. 6 (a) simplifies the concept of sensing using a satellite hyperspectral sensor. As can be seen, the sensor captures a pixel that contains three materials, i.e., Water, Tree, and Soil. At the sensor, irradiance is corrected to the reflectance in Fig. 6 (b). Note that this reflectance is between zero and one. Fig. 6 (c) shows the concept of linear unmixing, which demonstrates that the measured pixel contains 20% water (blue endmember), 50% soil (yellow endmember), and 30% tree (green endmember). Arguably, these portions cannot be negative, known as abundance nonnegativity constraint (ANC), and the summation of those portions is 100%, known as abundance sum to one constraint (ASC). Additionally, the endmembers are reflectance and cannot be negative (indeed, the endmembers have to be bounded between 0 and 1), which led to the endmember nonnegativity constraint. Here, we discuss two widely used linear mixing models, i.e., the low-rank linear mixing model (so-called LMM) and the sparse and redundant linear model. The former was mainly used for supervised and blind unmixing and the latter for sparse unmixing. Some other linear models have been used in the
<table><caption>TABLE I
THE LIST OF UNMIXING METHODS INCLUDED IN HYSUPP WITH THEIR CORRESPONDING DEPENDENCIES. WE HIGHLIGHTED OUR MAIN CONTRIBUTIONS AND THE LINK TO THE ORIGINAL IMPLEMENTATIONS.</caption><thead><tr><th>Method</th><th>Original implementation</th><th>Python</th><th>Dependencies</th><th>GPU</th><th>Contributions</th></tr></thead><tbody><tr><td>FCLSU [22]</td><td>pysptools</td><td>✓</td><td>numpy, cvxopt</td><td></td><td>Refactor into a SupervisedUnmixingModel sub-class</td></tr><tr><td>SiVM [23]</td><td>github</td><td>✓</td><td>numpy</td><td></td><td>Refactor into a BaseExtractor sub-class</td></tr><tr><td>SISAL [24]</td><td>github</td><td>✓</td><td>numpy</td><td></td><td>Refactor into a BaseExtractor sub-class</td></tr><tr><td>UnDIP [25]</td><td>github</td><td>✓</td><td>torch</td><td>✓</td><td>Refactor separate scripts into a single model</td></tr><tr><td>VCA [26]</td><td>github</td><td>✓</td><td> numpy</td><td></td><td>Replace scipy dependency by numpy. Add random seed</td></tr><tr><td>CLSLnSAL [27]</td><td>github</td><td>✓</td><td> numpy</td><td></td><td>Implementation based on SUNSAL</td></tr><tr><td>MUA_SLIC [28]</td><td>github</td><td>✓</td><td> numpy, skimage</td><td></td><td>MATLAB code translation using skimage's SLIC</td></tr><tr><td>S<sup>2</sup>WSU [29]</td><td>github</td><td>✓</td><td> numpy, scipy</td><td></td><td>MATLAB code translation</td></tr><tr><td>SUnAA [30]</td><td>github</td><td>✓</td><td> numpy, spams</td><td></td><td>Refactor to match SemiUnmixingModel sub-class</td></tr><tr><td>SUnCNN [21]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor separate scripts into a single model</td></tr><tr><td>SUnSAL [31]</td><td>github</td><td>✓</td><td> numpy</td><td></td><td>Refactor into SemiUnmixingModel sub-class</td></tr><tr><td>ADMMNet [32]</td><td>-</td><td>✓</td><td> torch</td><td>✓</td><td>Implemented from scratch</td></tr><tr><td>BayesianSMA [33]</td><td>webpage</td><td> matlab.engine</td><td> numpy</td><td></td><td>Python wrapper around existing MATLAB code</td></tr><tr><td>CNNAEU [34]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Convert existing keras implementation into torch</td></tr><tr><td>EDAA [35]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor to match BlindUnmixingModel sub-class</td></tr><tr><td>MiSiCNet [36]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor separate scripts into a single model</td></tr><tr><td>MSNet [37]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor to match BlindUnmixingModel sub-class</td></tr><tr><td>NMFQMV [38]</td><td>github</td><td> matlab.engine</td><td> numpy</td><td></td><td>Python wrapper around existing MATLAB code</td></tr><tr><td>PGMSU [39]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor to match BlindUnmixingModel sub-class</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Schematic of sensing and linear unmixing. (a) Sensing a mixed pixel (b) Reflectance of the mixed pixel and (c) Schematic of the linear unmixing.</figcaption></figure>

literature, which will be discussed throughout the paper. We also discuss the components of a basic neural network model.

Unmixing techniques can be categorized into three main groups considering the prior knowledge of endmembers. Supervised and unsupervised (blind) unmixing use the same low-rank mixing model, while semisupervised unmixing uses the sparse and redundant linear model. The endmember library is overcomplete; consequently, the abundances are desired to be sparse. Blind unmixing estimates both endmembers and abundances simultaneously. Fig. 7 graphically compares those three groups. It is worth mentioning that unsupervised (blind) unmixing can be applied to at sensor irradiance, however, the data cannot be interpreted and the endmembers extracted/estimated cannot be associated with the corresponding materials. One can apply supervised unmixing to irradiance if the pure pixels are also irradiance. Semisupervised unmixing cannot be applied to sensor radiance since they rely on a spectral library.

A. Low-rank Linear Mixture Model

Every $p$-dimensional pixel $\mathbf{y}$ (i.e., the sensor has $p$ bands) captured by a hyperspectral sensor can be represented as a linear combination of the endmembers within the pixel. Let us assume that matrix $\mathbf{E} \in \mathbb{R}^{p \times r}$ contains $r$ endmembers within the pixel. Then, $\mathbf{y} \in \mathbb{R}^p$ is represented as

$$ \mathbf{y} = \mathbf{E}\mathbf{a} + \mathbf{n}, \quad \text{s.t.} \quad \sum_{i=1}^{r} a_i = 1, a_i \ge 0, i = 1, 2, \dots, r, \quad (1) $$

where $\mathbf{n}$ denotes the $p$-dimensional random vector denoting the additive random Gaussian noise. Using matrix notations, we can represent all the pixels $\mathbf{Y}$ as

$$ \mathbf{Y} = \mathbf{E}\mathbf{A} + \mathbf{N}, \quad \text{s.t.} \quad \mathbf{A} \ge 0, 1_r^T \mathbf{A} = 1_n^T, \quad (2) $$

where $\mathbf{Y} \in \mathbb{R}^{p \times n}$ is the observed HSI, with $n$ pixels and $p$ bands, $\mathbf{N} \in \mathbb{R}^{p \times n}$ is noise, and $\mathbf{A} \in \mathbb{R}^{r \times n}$ contain the $r$ endmembers and their fractional abundances, respectively. $1_n$ indicates an $n$-component column vector of ones. Model 2 is known as the linear mixture model (LMM).
<table><caption>TABLE I
THE LIST OF UNMIXING METHODS INCLUDED IN HYSUPP WITH THEIR CORRESPONDING DEPENDENCIES. WE HIGHLIGHTED OUR MAIN CONTRIBUTIONS AND THE LINK TO THE ORIGINAL IMPLEMENTATIONS.</caption><thead><tr><th>Method</th><th>Original implementation</th><th>Python</th><th>Dependencies</th><th>GPU</th><th>Contributions</th></tr></thead><tbody><tr><td>FCLSU [22]</td><td>pysptools</td><td>✓</td><td>numpy, cvxopt</td><td></td><td>Refactor into a SupervisedUnmixingModel sub-class</td></tr><tr><td>SiVM [23]</td><td>github</td><td>✓</td><td>numpy</td><td></td><td>Refactor into a BaseExtractor sub-class</td></tr><tr><td>SISAL [24]</td><td>github</td><td>✓</td><td>numpy</td><td></td><td>Refactor into a BaseExtractor sub-class</td></tr><tr><td>UnDIP [25]</td><td>github</td><td>✓</td><td>torch</td><td>✓</td><td>Refactor separate scripts into a single model</td></tr><tr><td>VCA [26]</td><td>github</td><td>✓</td><td> numpy</td><td></td><td>Replace scipy dependency by numpy. Add random seed</td></tr><tr><td>CLSUnSAL [27]</td><td>github</td><td>✓</td><td> numpy</td><td></td><td>Implementation based on SUNSAL</td></tr><tr><td>MUA_SLAVE [28]</td><td>github</td><td>✓</td><td> numpy, skimage</td><td></td><td>MATLAB code translation using skimage's SLIC</td></tr><tr><td>S²WSU [29]</td><td>github</td><td>✓</td><td> numpy, scipy</td><td></td><td>MATLAB code translation</td></tr><tr><td>SUnAA [30]</td><td>github</td><td>✓</td><td> numpy, spams</td><td></td><td>Refactor to match SemiUnmixingModel sub-class</td></tr><tr><td>SUnCNN [21]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor separate scripts into a single model</td></tr><tr><td>SUnSAL [31]</td><td>github</td><td>✓</td><td> numpy</td><td></td><td>Refactor into SemiUnmixingModel sub-class</td></tr><tr><td>ADMMNet [32]</td><td>-</td><td>✓</td><td> torch</td><td>✓</td><td>Implemented from scratch</td></tr><tr><td>BayesianSMA [33]</td><td>webpage</td><td> matlab.engine</td><td> numpy</td><td></td><td>Python wrapper around existing MATLAB code</td></tr><tr><td>CNNAEU [34]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Convert existing keras implementation into torch</td></tr><tr><td>EDAA [35]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor to match BlindUnmixingModel sub-class</td></tr><tr><td>MiSiCNet [36]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor separate scripts into a single model</td></tr><tr><td>MSNet [37]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor to match BlindUnmixingModel sub-class</td></tr><tr><td>NMFQMV [38]</td><td>github</td><td> matlab.engine</td><td> numpy</td><td></td><td>Python wrapper around existing MATLAB code</td></tr><tr><td>PGMSU [39]</td><td>github</td><td>✓</td><td> torch</td><td>✓</td><td>Refactor to match BlindUnmixingModel sub-class</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Schematic of sensing and linear unmixing. (a) Sensing a mixed pixel (b) Reflectance of the mixed pixel and (c) Schematic of the linear unmixing.</figcaption></figure>

literature, which will be discussed throughout the paper. We also discuss the components of a basic neural network model.

Unmixing techniques can be categorized into three main groups considering the prior knowledge of endmembers. Supervised and unsupervised (blind) unmixing use the same low-rank mixing model, while semisupervised unmixing uses the sparse and redundant linear model. The endmember library is overcomplete; consequently, the abundances are desired to be sparse. Blind unmixing estimates both endmembers and abundances simultaneously. Fig. 7 graphically compares those three groups. It is worth mentioning that unsupervised (blind) unmixing can be applied to at sensor irradiance, however, the data cannot be interpreted and the endmembers extracted/estimated cannot be associated with the corresponding materials. One can apply supervised unmixing to irradiance if the pure pixels are also irradiance. Semisupervised unmixing cannot be applied to sensor radiance since they rely on a spectral library.

A. Low-rank Linear Mixture Model

Every p-dimensional pixel y (i.e., the sensor has p bands) captured by a hyperspectral sensor can be represented as a linear combination of the endmembers within the pixel. Let us assume that matrix E ∈ ℝ$^{p×r}$ contains r endmembers within the pixel. Then, y ∈ ℝ$^p$ is represented as

$$y = Ea + n, \quad \text{s.t.} \quad \sum_{i=1}^{r} a_i = 1, a_i \ge 0, i = 1, 2, \dots, r, \quad (1)$$

where n denotes the p-dimensional random vector denoting the additive random Gaussian noise. Using matrix notations, we can represent all the pixels Y as

$$Y = EA + N, \quad \text{s.t.} \quad A \ge 0, 1_r^T A = 1_n^T, \quad (2)$$

where Y ∈ ℝ$^{p×n}$ is the observed HSI, with n pixels and p bands, N ∈ ℝ$^{p×n}$ is noise, and A ∈ ℝ$^{r×n}$ contain the r endmembers and their fractional abundances, respectively. 1$_n$ indicates an n-component column vector of ones. Model 2 is known as the linear mixture model (LMM).
<figure><img src="image_8.png" /><figcaption>Fig. 7. Graphical representation of different types of linear unmixing. Supervised and blind unmixing use the same low-rank mixing model, while sparse unmixing uses the sparse and redundant linear model. Blind unmixing estimates both endmembers and abundances simultaneously.</figcaption></figure>

## B. Low-rank LMM and Simplex Volume

Suppose that the endmembers are affinely independent i.e., $\mathbf{e}_2 - \mathbf{e}_1, \dots, \mathbf{e}_r - \mathbf{e}_1$ are linearly independent, then

$$ S \triangleq \{\mathbf{E}\mathbf{a} \in \mathbb{R}^p | \sum_{i=1}^r a_i = 1, a_i \ge 0, i = 1, 2, \dots, r\} \quad (3) $$

is $(r-1)$-Simplex. Indeed, this is the convex hull of the vertices, i.e., $\mathbf{e}_i$'s assuming no noise. Therefore, if we ignore the noise, any point of the dataset belongs to $S$ (see Fig. 11 (a)). We should note that columns of $\mathbf{A}$ belong to $(r-1)$-probability simplex

$$ \Delta_r \triangleq \{\mathbf{a} \in \mathbb{R}^r | \sum_{i=1}^r a_i = 1, a_i \ge 0, i = 1, 2, \dots, r\}. \quad (4) $$

Therefore, $\mathbf{a}_i \in \Delta_r$ means ASC and ANC. Indeed, the vertices of $\Delta_r$ are $r$ unit vectors in $\mathbb{R}^r$. Hereafter, we refer to (3) as data simplex and to (4) as abundance simplex. It is worth mentioning that, unlike data simplex, the abundance simplex is also valid for nonlinear models.

Let us assume $S_r$ is an $r$-simplex in $\mathbb{R}^r$ i.e., $\mathbf{e}_i \in \mathbb{R}^r$ for $i = 0, 1, 2, \dots, r$. Then the volume of this simplex is given by

$$ \mathrm{Vol}(S_r) = \frac{1}{r!} |\det[\mathbf{e}_1 - \mathbf{e}_0, \dots, \mathbf{e}_r - \mathbf{e}_0]| = \frac{1}{r!} |\det \begin{bmatrix} 1 & \dots & 1 \\ \mathbf{e}_0 & \dots & \mathbf{e}_r \end{bmatrix}| \quad (5) $$

where $\det$ is the determinant of a matrix. In the linear mixture model (LMM), we have an $(r-1)$-Simplex in $\mathbb{R}^r$ (after the projection into an $r$ dimensional subspace, see IV-A) i.e., $\mathbf{e}_i \in \mathbb{R}^r$ for $i = 1, 2, \dots, r$ ($\mathbf{E}_r \in \mathbb{R}^{r \times r}$). The volume of an $(r-1)$-simplex is zero in $\mathbb{R}^r$. Therefore, the extended simplex, which contains the origin (i.e., $\mathbf{E}_0 = [0, \mathbf{E}]$), can be used to calculate the volume. Hence,

$$ \mathrm{Vol}(\mathbf{E}_0) = \frac{1}{r!} |\det[\mathbf{e}_1, \dots, \mathbf{e}_r]| = \frac{1}{r!} |\det(\mathbf{E})|. \quad (6) $$

Alternatively, the data points can be shifted to the origin, and the volume of the $(r-1)$-simplex can be computed in $\mathbb{R}^{r-1}$ given by

$$ \mathrm{Vol}(\mathbf{E}_r) = \frac{1}{(r-1)!} |\det[\mathbf{e}_2 - \mathbf{e}_1, \dots, \mathbf{e}_r - \mathbf{e}_1]| = \frac{1}{(r-1)!} |\det \begin{bmatrix} 1 & \dots & 1 \\ \mathbf{e}_1 & \dots & \mathbf{e}_r \end{bmatrix}|. \quad (7) $$

It is worth mentioning that the projection into a subspace (IV-A) is necessary to form a squared matrix and calculate the determinant.

## C. Sparse and Redundant Linear Mixture Model

The observed spectra can be represented using a sparse and redundant linear mixture model given by

$$ \mathbf{Y} = \mathbf{D}\mathbf{X} + \mathbf{N}, \\ \text{s.t.} \quad \mathbf{X} \ge 0, \mathbf{1}_m^T\mathbf{X} = \mathbf{1}_n^T, \quad (8) $$

where $\mathbf{D} \in \mathbb{R}^{p \times m}$ ($p \ll m$) denotes the spectral library containing $m$ endmembers and $\mathbf{X} \in \mathbb{R}^{m \times n}$ is the unknown fractional abundances to estimate. Note that $\mathbf{D}$ is an over-complete dictionary and therefore, should be a well-designed dictionary. A well-designed dictionary contains endmembers of the materials in the scene and can sparsify the redundant $\mathbf{X}$. Hence, a spectral library can be pruned based on the spectral angles of the spectra (i.e., spectra with small degree differences will be removed). However, this comes with a risk of losing material endmembers if they are scaled versions of each other. In the case of a well-designed dictionary, the pixels are a mixture of a few atoms of the dictionary, and therefore $\mathbf{X}$ is a sparse matrix. Note that we get a zero row if an endmember material does not exist in the observed spectra. The rest of the rows will contain zeros since the abundances are often sparse. This model is often used in sparse unmixing. The fractional abundances $\mathbf{X}$ are estimated using sparsity-enforcing penalties/constraints in a sparse regression formulation.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Graphical representation of different types of linear unmixing. Supervised and blind unmixing use the same low-rank mixing model, while sparse unmixing uses the sparse and redundant linear model. Blind unmixing estimates both endmembers and abundances simultaneously.</figcaption></figure>

## B. Low-rank LMM and Simplex Volume

Suppose that the endmembers are affinely independent i.e., $\mathbf{e}_2 - \mathbf{e}_1, ..., \mathbf{e}_r - \mathbf{e}_1$ are linearly independent, then

$$ S \triangleq \{\mathbf{E}\mathbf{a} \in \mathbb{R}^p | \sum_{i=1}^r a_i = 1, a_i \ge 0, i = 1, 2, \dots, r\} \quad (3) $$

is $(r-1)$-Simplex. Indeed, this is the convex hull of the vertices, i.e., $\mathbf{e}_i$'s assuming no noise. Therefore, if we ignore the noise, any point of the dataset belongs to $S$ (see Fig. 11 (a)). We should note that columns of $\mathbf{A}$ belong to $(r-1)$-probability simplex

$$ \Delta_r \triangleq \{\mathbf{a} \in \mathbb{R}^r | \sum_{i=1}^r a_i = 1, a_i \ge 0, i = 1, 2, \dots, r\}. \quad (4) $$

Therefore, $\mathbf{a}_i \in \Delta_r$ means ASC and ANC. Indeed, the vertices of $\Delta_r$ are $r$ unit vectors in $\mathbb{R}^r$. Hereafter, we refer to (3) as data simplex and to (4) as abundance simplex. It is worth mentioning that, unlike data simplex, the abundance simplex is also valid for nonlinear models.

Let us assume $S_r$ is an $r$-simplex in $\mathbb{R}^r$ i.e., $\mathbf{e}_i \in \mathbb{R}^r$ for $i = 0, 1, 2, \dots, r$. Then the volume of this simplex is given by

$$ \mathrm{Vol}(S_r) = \frac{1}{r!} |\det[\mathbf{e}_1 - \mathbf{e}_0, \dots, \mathbf{e}_r - \mathbf{e}_0]| = \frac{1}{r!} |\det \begin{bmatrix} 1 & \dots & 1 \\ \mathbf{e}_0 & \dots & \mathbf{e}_r \end{bmatrix}| \quad (5) $$

where $\det$ is the determinant of a matrix. In the linear mixture model (LMM), we have an $(r-1)$-Simplex in $\mathbb{R}^r$ (after the projection into an $r$ dimensional subspace, see IV-A) i.e., $\mathbf{e}_i \in \mathbb{R}^r$ for $i = 1, 2, \dots, r$ ($\mathbf{E}_r \in \mathbb{R}^{r \times r}$). The volume of an $(r-1)$-simplex is zero in $\mathbb{R}^r$. Therefore, the extended simplex, which contains the origin (i.e., $\mathbf{E}_0 = [0, \mathbf{E}]$), can be used to calculate the volume. Hence,

$$ \mathrm{Vol}(\mathbf{E}_0) = \frac{1}{r!} |\det[\mathbf{e}_1, \dots, \mathbf{e}_r]| = \frac{1}{r!} |\det(\mathbf{E})|. \quad (6) $$

Alternatively, the data points can be shifted to the origin, and the volume of the $(r-1)$-simplex can be computed in $\mathbb{R}^{r-1}$ given by

$$ \mathrm{Vol}(\mathbf{E}_r) = \frac{1}{(r-1)!} |\det[\mathbf{e}_2 - \mathbf{e}_1, \dots, \mathbf{e}_r - \mathbf{e}_1]| = \frac{1}{(r-1)!} |\det \begin{bmatrix} 1 & \dots & 1 \\ \mathbf{e}_1 & \dots & \mathbf{e}_r \end{bmatrix}|. \quad (7) $$

It is worth mentioning that the projection into a subspace (IV-A) is necessary to form a squared matrix and calculate the determinant.

## C. Sparse and Redundant Linear Mixture Model

The observed spectra can be represented using a sparse and redundant linear mixture model given by

$$ \mathbf{Y} = \mathbf{D}\mathbf{X} + \mathbf{N}, \\ \text{s.t.} \quad \mathbf{X} \ge 0, \mathbf{1}_m^T\mathbf{X} = \mathbf{1}_n^T, \quad (8) $$

where $\mathbf{D} \in \mathbb{R}^{p \times m}$ ($p \ll m$) denotes the spectral library containing $m$ endmembers and $\mathbf{X} \in \mathbb{R}^{m \times n}$ is the unknown fractional abundances to estimate. Note that $\mathbf{D}$ is an over-complete dictionary and therefore, should be a well-designed dictionary. A well-designed dictionary contains endmembers of the materials in the scene and can sparsify the redundant $\mathbf{X}$. Hence, a spectral library can be pruned based on the spectral angles of the spectra (i.e., spectra with small degree differences will be removed). However, this comes with a risk of losing material endmembers if they are scaled versions of each other. In the case of a well-designed dictionary, the pixels are a mixture of a few atoms of the dictionary, and therefore $\mathbf{X}$ is a sparse matrix. Note that we get a zero row if an endmember material does not exist in the observed spectra. The rest of the rows will contain zeros since the abundances are often sparse. This model is often used in sparse unmixing. The fractional abundances $\mathbf{X}$ are estimated using sparsity-enforcing penalties/constraints in a sparse regression formulation.
## D. Shallow/Deep Neural Network Model

Shallow or deep neural networks can be used for supervised, semi-supervised, or blind unmixing. The network is often based on autoencoder-based architecture, as shown in Fig. 8. The encoder $\mathcal{E}$, which could be deep or shallow, encodes the spectral pixels into the abundances given by

$$ \mathbf{a} = \mathcal{E}(\mathbf{y}), \qquad (9) $$

and the decoder $\mathcal{D}$ reconstruct the pixel given by

$$ \hat{\mathbf{y}} = \mathcal{D}(\mathbf{a}). \qquad (10) $$

In every layer, an activation function is used to apply non-linearity. The most common choices are Rectified Linear Unit ReLU (ReLU($x$) = $max(0, x)$) and Leaky ReLU (LReLU($x$) = $max(ax, x)$, where $a$ is very small, and a common choice is $a = 0.1$). Batch normalization (BN) is commonly used in each layer to speed up the learning process. The BN function is given by

$$ BN(x) = \alpha_2 \frac{x - \mu}{\sqrt{\sigma^2 + \epsilon}} + \alpha_1, \qquad (11) $$

where $\alpha_2$ and $\alpha_1$ are learnable parameters, $\mu$ and $\sigma$ are the mean and standard deviation of the batch, and $\epsilon$ is a minimal value.

A linear layer can be used for the decoder to reconstruct the pixel. The decoder weights are the endmembers fixed for supervised unmixing and can be learned through the learning process in the case of unsupervised (blind) unmixing. ASC and ANC can be realized differently; Applying nonnegative activation functions such as ReLU and normalization for ASC given by [10]

$$ \mathbf{a} = \frac{\mathbf{a}}{\sum_{i=1}^{r} a_i}, \qquad (12) $$

alternatively, one can use

$$ \mathbf{a} = \frac{|\mathbf{a}|}{\sum_{i=1}^{r} |a_i|}, \qquad (13) $$

A common way is to use a softmax function

$$ \mathrm{softmax}(\mathbf{a}) = \frac{e^{\mathbf{a}}}{\sum_{i=1}^{r} e^{a_i}}. \qquad (14) $$

The scaled version of the softmax can be used as $softmax(\gamma\mathbf{a}_i)$ where $\gamma > 1$ is constant. ASC and ANC can also be implemented by adding a penalty term to the loss function. However, it degrades the stability of the training process and decreases the convergence time, and therefore, it is not recommended [10]. We should note that those functions are realized as the last hidden layer of the decoder. The loss function is often chosen as a reconstruction term plus a regularization term given by

$$ \mathcal{L}_{\mathrm{reg}}(\mathbf{y}, \hat{\mathbf{y}}) = \mathcal{L}(\mathbf{y}, \hat{\mathbf{y}}) + \lambda J(\mathbf{a}, \hat{\mathbf{W}}_{\mathrm{en}}, \hat{\mathbf{W}}_{\mathrm{de}}), \qquad (15) $$

where $\mathcal{L}$ is often selected as a reconstruction term such as $\ell_2$ norm (or MSE)

$$ \mathcal{L}_{\ell_2}(\mathbf{y}, \hat{\mathbf{y}}) = \|\mathbf{y} - \hat{\mathbf{y}}\|_2^2, \qquad (16) $$

or

$$ \mathcal{L}_{\mathrm{SAD}}(\mathbf{y}, \hat{\mathbf{y}}) = \arccos \frac{\mathbf{y}^T \hat{\mathbf{y}}}{\|\mathbf{y}\|_2 \|\hat{\mathbf{y}}\|_2}, \qquad (17) $$

or

$$ \mathcal{L}_{\mathrm{SID}}(\mathbf{y}, \hat{\mathbf{y}}) = \sum_{i=1}^{p} p_i \log\left(\frac{p_i}{q_i}\right) + \sum_{i=1}^{p} q_i \log\left(\frac{q_i}{p_i}\right), \qquad (18) $$

where

$$ p_i = \frac{y_i}{\sum_{i=1}^{p} y_i}, \quad q_i = \frac{\hat{y}_i}{\sum_{i=1}^{p} \hat{y}_i}, \qquad (19) $$

or a combination of them. The regularizer $J$ can be a combination of several functions, such as a Frobenius norm on $\mathbf{W}_{\mathrm{en}}$ with a sparsity norm such as $\ell_1$ norm on $\mathbf{a}$, and a Frobenius norm or a geometrical (MV induced) norm on $\mathbf{W}_{\mathrm{de}}$. Note that SAD is scale-invariant and, therefore, can capture the spectral variability of the pixel. This is a big advantage for real-world datasets. On the other hand, utilizing SAD for datasets containing materials having scaled endmembers is not recommended since the network cannot be trained to distinguish those materials. Assuming a single-layer encoder-decoder, the abundances are given by

$$ \mathbf{a} = \mathrm{softmax}(BN(\mathrm{ReLU}(\mathbf{W}_{\mathrm{en}}\mathbf{y} + \mathbf{b}))), \qquad (20) $$

where $b$ denotes the biases of the layer. The reconstructed pixel is given by

$$ \hat{\mathbf{y}} = \mathbf{W}_{\mathrm{de}}\mathbf{a}. \qquad (21) $$

<figure><img src="image_9.png" /><figcaption>Fig. 8. A basic encoder-decode architecture used for unmixing. The (shallow or deep) encoder encodes the pixels to abundance. ASC and ANC can be enforced on the bottleneck. The shallow decoder reconstructs the pixels, and the decoder's weights are the endmembers.</figcaption></figure>

## E. The Effect of Noise in Unmixing

Hyperspectral data are usually downgraded by different sources, including atmospheric effects, imaging artifacts, and instrumental noise. These sources can further affect hyperspectral analysis, including unmixing [5]. Unmixing is an inverse problem that could be very sensitive to noise, mainly when endmembers are highly correlated. Therefore, unmixing techniques may fail due to noise. In [1], so-called signal-to-noise-ratio spectral distribution (SNR-SD) was suggested to determine if the unmixing inverse problem gives acceptable results. However, the hyperspectral denoising field has considerably evolved during the past decade [4], [5], [40] and noise reduction as a preprocessing step could improve the unmixing performance [41]. We should note that unmixing techniques usually consider Gaussian noise and even can perform as a denoiser but are not as efficient as denoising techniques [41]. Additionally, other types of noise, such as
## D. Shallow/Deep Neural Network Model

Shallow or deep neural networks can be used for supervised, semi-supervised, or blind unmixing. The network is often based on autoencoder-based architecture, as shown in Fig. 8. The encoder $\mathcal{E}$, which could be deep or shallow, encodes the spectral pixels into the abundances given by

$$ \mathbf{a} = \mathcal{E}(\mathbf{y}), \qquad (9) $$

and the decoder $\mathcal{D}$ reconstruct the pixel given by

$$ \hat{\mathbf{y}} = \mathcal{D}(\mathbf{a}). \qquad (10) $$

In every layer, an activation function is used to apply non-linearity. The most common choices are Rectified Linear Unit ReLU (ReLU($x$) = $max(0, x)$) and Leaky ReLU (LReLU($x$) = $max(ax, x)$), where $a$ is very small, and a common choice is $a = 0.1$. Batch normalization (BN) is commonly used in each layer to speed up the learning process. The BN function is given by

$$ BN(x) = \alpha_2 \frac{x - \mu}{\sqrt{\sigma^2 + \epsilon}} + \alpha_1, \qquad (11) $$

where $\alpha_2$ and $\alpha_1$ are learnable parameters, $\mu$ and $\sigma$ are the mean and standard deviation of the batch, and $\epsilon$ is a minimal value.

A linear layer can be used for the decoder to reconstruct the pixel. The decoder weights are the endmembers fixed for supervised unmixing and can be learned through the learning process in the case of unsupervised (blind) unmixing. ASC and ANC can be realized differently; Applying nonnegative activation functions such as ReLU and normalization for ASC given by [10]

$$ \mathbf{a} = \frac{\mathbf{a}}{\sum_{i=1}^{r} a_i}, \qquad (12) $$

alternatively, one can use

$$ \mathbf{a} = \frac{|\mathbf{a}|}{\sum_{i=1}^{r} |a_i|}, \qquad (13) $$

A common way is to use a softmax function

$$ \mathrm{softmax}(\mathbf{a}) = \frac{e^{\mathbf{a}}}{\sum_{i=1}^{r} e^{a_i}}. \qquad (14) $$

The scaled version of the softmax can be used as $softmax(\gamma\mathbf{a}_i)$ where $\gamma > 1$ is constant. ASC and ANC can also be implemented by adding a penalty term to the loss function. However, it degrades the stability of the training process and decreases the convergence time, and therefore, it is not recommended [10]. We should note that those functions are realized as the last hidden layer of the decoder. The loss function is often chosen as a reconstruction term plus a regularization term given by

$$ \mathcal{L}_{reg}(\mathbf{y}, \hat{\mathbf{y}}) = \mathcal{L}(\mathbf{y}, \hat{\mathbf{y}}) + \lambda J(\mathbf{a}, \hat{\mathbf{W}}_{en}, \hat{\mathbf{W}}_{de}), \qquad (15) $$

where $\mathcal{L}$ is often selected as a reconstruction term such as $\ell_2$ norm (or MSE)

$$ \mathcal{L}_{\ell_2}(\mathbf{y}, \hat{\mathbf{y}}) = \|\mathbf{y} - \hat{\mathbf{y}}\|_2^2, \qquad (16) $$

or

$$ \mathcal{L}_{SAD}(\mathbf{y}, \hat{\mathbf{y}}) = \arccos \frac{\mathbf{y}^T \hat{\mathbf{y}}}{\|\mathbf{y}\|_2 \|\hat{\mathbf{y}}\|_2}, \qquad (17) $$

or

$$ \mathcal{L}_{SID}(\mathbf{y}, \hat{\mathbf{y}}) = \sum_{i=1}^{p} p_i \log\left(\frac{p_i}{q_i}\right) + \sum_{i=1}^{p} q_i \log\left(\frac{q_i}{p_i}\right), \qquad (18) $$

where

$$ p_i = \frac{y_i}{\sum_{i=1}^{p} y_i}, \quad q_i = \frac{\hat{y}_i}{\sum_{i=1}^{p} \hat{y}_i}, \qquad (19) $$

or a combination of them. The regularizer $J$ can be a combination of several functions, such as a Frobenius norm on $\mathbf{W}_{en}$ with a sparsity norm such as $\ell_1$ norm on $\mathbf{a}$, and a Frobenius norm or a geometrical (MV induced) norm on $\mathbf{W}_{de}$. Note that SAD is scale-invariant and, therefore, can capture the spectral variability of the pixel. This is a big advantage for real-world datasets. On the other hand, utilizing SAD for datasets containing materials having scaled endmembers is not recommended since the network cannot be trained to distinguish those materials. Assuming a single-layer encoder-decoder, the abundances are given by

$$ \mathbf{a} = \mathrm{softmax}(BN(\mathrm{RelU}(\mathbf{W}_{en}\mathbf{y} + \mathbf{b}))), \qquad (20) $$

where $b$ denotes the biases of the layer. The reconstructed pixel is given by

$$ \hat{\mathbf{y}} = \mathbf{W}_{de}\mathbf{a}. \qquad (21) $$

<figure><img src="image_9.png" /><figcaption>Fig. 8. A basic encoder-decode architecture used for unmixing. The (shallow or deep) encoder encodes the pixels to abundance. ASC and ANC can be enforced on the bottleneck. The shallow decoder reconstructs the pixels, and the decoder's weights are the endmembers.</figcaption></figure>

## E. The Effect of Noise in Unmixing

Hyperspectral data are usually downgraded by different sources, including atmospheric effects, imaging artifacts, and instrumental noise. These sources can further affect hyperspectral analysis, including unmixing [5]. Unmixing is an inverse problem that could be very sensitive to noise, mainly when endmembers are highly correlated. Therefore, unmixing techniques may fail due to noise. In [1], so-called signal-to-noise-ratio spectral distribution (SNR-SD) was suggested to determine if the unmixing inverse problem gives acceptable results. However, the hyperspectral denoising field has considerably evolved during the past decade [4], [5], [40] and noise reduction as a preprocessing step could improve the unmixing performance [41]. We should note that unmixing techniques usually consider Gaussian noise and even can perform as a denoiser but are not as efficient as denoising techniques [41]. Additionally, other types of noise, such as
<figure><img src="image_10.png" /><figcaption>Fig. 9. Taxonomy of the linear unmixing methods described in this paper. The top categorization is based on the prior knowledge of endmembers.</figcaption></figure>

strip and sparse, can be removed by applying denoising before unmixing [40]. Alternatively, some methods were proposed for performing denoising and unmixing in a unified framework for boosting the performance of each other [42], [43]. As a result, we suggest applying a well-established noise reduction technique before estimating the abundance or characterizing the unmixing inverse problem using SNR-SD.

## F. Unmixing Approaches w.r.t. Prior Knowledge of Endmembers

Unmixing techniques can be categorized in terms of prior knowledge of the endmembers. Depending on the prior knowledge available about endmembers, the unmixing problem can be divided into three main categories: (1) supervised unmixing, (2) unsupervised (blind) unmixing, and (3) Semi-supervised unmixing. In supervised unmixing, abundances are estimated by relying on known endmembers, whereas blind unmixing estimates both the endmembers and the abundances simultaneously. Semi-supervised unmixing relies on a library of endmembers that ideally contains the endmembers in the scene. We discuss each category separately. Fig. 7 demonstrates the graphical abstracts of these three linear unmixing categories. The taxonomy of the linear unmixing methods discussed in this paper is given in Fig. 9.

# IV. SUPERVISED UNMIXING

In supervised unmixing, we assume that endmembers are known, and the abundance matrix needs to be estimated. In practice, the endmember can be measured in the field or laboratory. One can select them from a spectral library, however, due to the variations in the imaging setup, such a selection will not lead to a desirable abundance estimation. Alternatively,

they can be extracted/estimated from data points. Geometrical approaches are often used at this step. Endmember extraction is not often easy since the captured data might not contain pure pixels to represent all the endmembers in the scene. Most endmember extraction methods rely on pure pixels or some pixels on the facet of the simplex. Geometrical approaches are less effective in no pure pixel scenarios. The processing chain of supervised unmixing is shown in Fig. 10. Note that we classify methods based on such sequential process (as opposed to simultaneous estimation of endmembers and abundances) as supervised unmixing since in the step of abundance estimation the endmembers are assumed to be known. The processing chain often includes three main steps; 1) Subspace projection: the data is projected into a subspace, 2) endmember extraction/estimation: Often a geometrical approach is used to extract the endmembers, and 3) Abundance estimation: abundances can be estimated using a least square (LS) or a neural network (NN) -based approach.

## A. Subspace Projection and Estimation of number of Endmembers

Hyperspectral data often live in subspace having dimension much lower than the dimension of spectral bands defined by the sensor. Assuming *r* endmembers in the scene, the intrinsic/subspace dimension is *r* − 1, i.e., the data points can be represented by *r* − 1 linearly independent vectors or bases (in the case of orthogonal projections). Therefore, identifying such a subspace and projecting the data into that reduces the computational cost, memory consumption, and removes the noise and outliers.

The subspace projection is given by

$$ \mathbf{Y} = \mathbf{V}^{T}\mathbf{F} + \mathbf{N}, \qquad (22) $$
<figure><img src="image_10.png" /><figcaption>Fig. 9. Taxonomy of the linear unmixing methods described in this paper. The top categorization is based on the prior knowledge of endmembers.</figcaption></figure>

strip and sparse, can be removed by applying denoising before unmixing [40]. Alternatively, some methods were proposed for performing denoising and unmixing in a unified framework for boosting the performance of each other [42], [43]. As a result, we suggest applying a well-established noise reduction technique before estimating the abundance or characterizing the unmixing inverse problem using SNR-SD.

## F. Unmixing Approaches w.r.t. Prior Knowledge of Endmembers

Unmixing techniques can be categorized in terms of prior knowledge of the endmembers. Depending on the prior knowledge available about endmembers, the unmixing problem can be divided into three main categories: (1) supervised unmixing, (2) unsupervised (blind) unmixing, and (3) Semi-supervised unmixing. In supervised unmixing, abundances are estimated by relying on known endmembers, whereas blind unmixing estimates both the endmembers and the abundances simultaneously. Semi-supervised unmixing relies on a library of endmembers that ideally contains the endmembers in the scene. We discuss each category separately. Fig. 7 demonstrates the graphical abstracts of these three linear unmixing categories. The taxonomy of the linear unmixing methods discussed in this paper is given in Fig. 9.

# IV. SUPERVISED UNMIXING

In supervised unmixing, we assume that endmembers are known, and the abundance matrix needs to be estimated. In practice, the endmember can be measured in the field or laboratory. One can select them from a spectral library, however, due to the variations in the imaging setup, such a selection will not lead to a desirable abundance estimation. Alternatively,

they can be extracted/estimated from data points. Geometrical approaches are often used at this step. Endmember extraction is not often easy since the captured data might not contain pure pixels to represent all the endmembers in the scene. Most endmember extraction methods rely on pure pixels or some pixels on the facet of the simplex. Geometrical approaches are less effective in no pure pixel scenarios. The processing chain of supervised unmixing is shown in Fig. 10. Note that we classify methods based on such sequential process (as opposed to simultaneous estimation of endmembers and abundances) as supervised unmixing since in the step of abundance estimation the endmembers are assumed to be known. The processing chain often includes three main steps; 1) Subspace projection: the data is projected into a subspace, 2) endmember extraction/estimation: Often a geometrical approach is used to extract the endmembers, and 3) Abundance estimation: abundances can be estimated using a least square (LS) or a neural network (NN) -based approach.

## A. Subspace Projection and Estimation of number of Endmembers

Hyperspectral data often live in subspace having dimension much lower than the dimension of spectral bands defined by the sensor. Assuming *r* endmembers in the scene, the intrinsic/subspace dimension is *r* − 1, i.e., the data points can be represented by *r* − 1 linearly independent vectors or bases (in the case of orthogonal projections). Therefore, identifying such a subspace and projecting the data into that reduces the computational cost, memory consumption, and removes the noise and outliers.

The subspace projection is given by

$$ \mathbf{Y} = \mathbf{V}^T \mathbf{F} + \mathbf{N}, \qquad (22) $$
<figure><img src="image_11.png" /><figcaption>Fig. 10. The processing chain of supervised unmixing. First, the data is projected into a subspace. Then, endmembers are extracted using a geometrical approach. If the endmembers are available, then the endmember extraction can be skipped. In the final step, abundances can be estimated using an ML/DL-based approach.</figcaption></figure>

where columns $V \in \mathbb{R}^{p \times r}$ spans the subspace and columns of $F$ are the projected spectral pixels. When $V$ is semi orthogonal, i.e., $V^T V = I_r$, the project is orthogonal which is a common choice in HS subspace projection. In (22), both $V$ and $F$ should be estimated [44], [45]. Alternately, we can assume [1]

$$Y_r = V^T E A + V^T N, \quad (23)$$

where $Y_r = V^T Y$. Reduction methods such as minimum noise fraction (MNF), principal component analysis (PCA), and noise adjusted principal component (NAPC) estimate $V$ [46]. Not that $Y$ is usually not a square matrix; hence, singular value decomposition (SVD) can be used to estimate $V$ [47], [48]. Please see [16], [46] for more details on hyperspectral feature reduction.

Estimating the number of endmembers is not a trivial task. Unmixing approaches are vulnerable to this parameter and under or overestimating $r$ may considerably affect the error of the models. We should note that sparse regression-based methods which do not use endmember bundles do not rely on the estimation of $r$. In the literature, this problem was addressed with alternative names such as hyperspectral subspace identification, intrinsic order selection, virtual dimension, and estimation of the number of spectrally distinct signal sources [47]–[49]. This problem can be addressed using eigenvalue-based detection techniques [49]–[51] or estimating the mean square errors [47], [48]. Geometrical-based approaches were also proposed for endmember estimation [52]–[54].

## B. Endmember Extraction for Linear Unmixing

As can be seen in Fig. 11 (a), the endmembers are the vertices of the simplex enclosing the data points. Therefore, the geometrical concept inspires the mainstream of endmember extraction techniques. Geometrical approaches can be divided into two main groups; pure pixel-based approaches and approaches without pure pixel assumption. Fig. 11 compares the pure pixel with no pure pixel assumption. In Fig. 11 (b), the red circles demonstrate the positions of the missing pure pixels. In this case, the endmembers cannot be extracted and they should be estimated. This is feasible using geometrical approaches if enough pixels are located on the facets of the simplex. This becomes much more challenging for higher dimensional cases when noise and other nonlinearities are involved. Alternatively, a group of endmembers (endmember bundle) can be extracted based on the assumption of endmember variability caused by spectral variability as shown in Fig. 11 (c). Fig. 11 (d) shows

the highly mixed scenario that geometrical approaches cannot cope with.

<figure><img src="image_12.png" /><figcaption>Fig. 11. Pure pixel versus no pure pixel assumption. (a) Pure pixels scenario (b) No pure pixels scenario with enough pixels on the facets of the simplex (c) endmember bundle, and (d) No pure pixels scenario without pixels on the facets of the simplex (highly mixed scenario).</figcaption></figure>

1) Pure Pixel Assumption: Many techniques consider pure pixel assumption for simplicity, while in real applications, pure pixels of some endmembers are often missing. The methods that rely on pure pixels for endmember extraction can be divided into three main groups: projections and extremes, simplex fitting methods, and multiple endmember extraction methods (endmember bundles).

* Projections and Extremes: This group often searches for extremes by iteratively projecting data points. The vertices can be selected as the extreme points after iteratively projecting the data in a particular direction. For instance, Pixel Purity Index (PPI) [55] scores the spectral vectors by projecting them onto a large set of random vectors (called skewers) and counting the number of times that each vector is an extreme point. Orthogonal subspace projection (OSP) [56] and Vertex Component Analysis (VCA) [26] selects endmembers iteratively by projecting the data into an orthogonal direction to the subspace spanned by the already selected endmember.

* Simplex Fitting: The endmembers are assumed to be located at the vertices of the simplex enclosing the data points. Therefore, they can be extracted by maximizing the data simplex. N-FINDR [57] searches for pure pixels that form the largest simplex by gradually inflating a
<figure><img src="image_11.png" /><figcaption>Fig. 10. The processing chain of supervised unmixing. First, the data is projected into a subspace. Then, endmembers are extracted using a geometrical approach. If the endmembers are available, then the endmember extraction can be skipped. In the final step, abundances can be estimated using an ML/DL-based approach.</figcaption></figure>

where columns $V \in \mathbb{R}^{p \times r}$ spans the subspace and columns of $F$ are the projected spectral pixels. When $V$ is semi orthogonal, i.e., $V^T V = I_r$, the project is orthogonal which is a common choice in HS subspace projection. In (22), both $V$ and $F$ should be estimated [44], [45]. Alternately, we can assume [1]

$$Y_r = V^T E A + V^T N, \quad (23)$$

where $Y_r = V^T Y$. Reduction methods such as minimum noise fraction (MNF), principal component analysis (PCA), and noise adjusted principal component (NAPC) estimate $V$ [46]. Not that $Y$ is usually not a square matrix; hence, singular value decomposition (SVD) can be used to estimate $V$ [47], [48]. Please see [16], [46] for more details on hyperspectral feature reduction.

Estimating the number of endmembers is not a trivial task. Unmixing approaches are vulnerable to this parameter and under or overestimating $r$ may considerably affect the error of the models. We should note that sparse regression-based methods which do not use endmember bundles do not rely on the estimation of $r$. In the literature, this problem was addressed with alternative names such as hyperspectral subspace identification, intrinsic order selection, virtual dimension, and estimation of the number of spectrally distinct signal sources [47]–[49]. This problem can be addressed using eigenvalue-based detection techniques [49]–[51] or estimating the mean square errors [47], [48]. Geometrical-based approaches were also proposed for endmember estimation [52]–[54].

## B. Endmember Extraction for Linear Unmixing

As can be seen in Fig. 11 (a), the endmembers are the vertices of the simplex enclosing the data points. Therefore, the geometrical concept inspires the mainstream of endmember extraction techniques. Geometrical approaches can be divided into two main groups; pure pixel-based approaches and approaches without pure pixel assumption. Fig. 11 compares the pure pixel with no pure pixel assumption. In Fig. 11 (b), the red circles demonstrate the positions of the missing pure pixels. In this case, the endmembers cannot be extracted and they should be estimated. This is feasible using geometrical approaches if enough pixels are located on the facets of the simplex. This becomes much more challenging for higher dimensional cases when noise and other nonlinearities are involved. Alternatively, a group of endmembers (endmember bundle) can be extracted based on the assumption of endmember variability caused by spectral variability as shown in Fig. 11 (c). Fig. 11 (d) shows

the highly mixed scenario that geometrical approaches cannot cope with.

<figure><img src="image_12.png" /><figcaption>Fig. 11. Pure pixel versus no pure pixel assumption. (a) Pure pixels scenario (b) No pure pixels scenario with enough pixels on the facets of the simplex (c) endmember bundle, and (d) No pure pixels scenario without pixels on the facets of the simplex (highly mixed scenario).</figcaption></figure>

1) Pure Pixel Assumption: Many techniques consider pure pixel assumption for simplicity, while in real applications, pure pixels of some endmembers are often missing. The methods that rely on pure pixels for endmember extraction can be divided into three main groups: projections and extremes, simplex fitting methods, and multiple endmember extraction methods (endmember bundles).

* Projections and Extremes: This group often searches for extremes by iteratively projecting data points. The vertices can be selected as the extreme points after iteratively projecting the data in a particular direction. For instance, Pixel Purity Index (PPI) [55] scores the spectral vectors by projecting them onto a large set of random vectors (called skewers) and counting the number of times that each vector is an extreme point. Orthogonal subspace projection (OSP) [56] and Vertex Component Analysis (VCA) [26] selects endmembers iteratively by projecting the data into an orthogonal direction to the subspace spanned by the already selected endmember.

* Simplex Fitting: The endmembers are assumed to be located at the vertices of the simplex enclosing the data points. Therefore, they can be extracted by maximizing the data simplex. N-FINDR [57] searches for pure pixels that form the largest simplex by gradually inflating a
simplex inside the data. Simplex volume maximization (SiVM) [23] extracts the endmembers by iteratively maximizing the simplex volume using

$$ \arg\max_{\mathbf{E}_r} V(\mathbf{E}_r) = \arg\max_{\mathbf{E}_r} \sqrt{\frac{(-\mathbf{1})^r \cdot \mathrm{cmd}(\mathbf{E}_r)}{2^{r-1}(\mathbf{r}-1)!}}, \quad (24) $$

where cmd is the Cayley–Menger determinant:

$$ \mathrm{cmd}(\mathbf{E}_r) = \det \begin{bmatrix} 0 & 1 & 1 & 1 & \cdots & 1 \\ 1 & 0 & d_{1,2}^2 & d_{1,3}^2 & \cdots & d_{1,r}^2 \\ 1 & d_{2,1}^2 & 0 & d_{2,3}^2 & \cdots & d_{2,r}^2 \\ 1 & d_{3,1}^2 & d_{3,2}^2 & 0 & \cdots & d_{3,r}^2 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & d_{r,1}^2 & d_{r,2}^2 & d_{r,3}^2 & \cdots & 0 \end{bmatrix}, $$

and $d_{i,j}$ is the Euclidean distance between endmembers $e_i$ and $e_j$.

*   **Multiple Pure Pixels and Endmember Bundles:** The multiple endmembers were also referred to as the endmember bundles. The endmember bundle and bundle unmixing were proposed in [58]–[60]. PCA was used to identify the (projected) endmembers using the extremes. More than one extreme was selected for materials. The bundle unmixing was formed using linear programming to determine the minimum, average, and maximum of the fractional abundances induced by endmember variability. In [61], endmember bundles were extracted from the observed data. Endmember variability is modeled using the endmember bundle, which will be discussed in more detail in semisupervised unmixing.

Some techniques incorporate spatial information for extracting the endmembers [62]–[65]. These methods often utilize information from neighboring pixels in the endmember extraction. For instance, in [66], the similarity of the pixel and its neighbors was used to weigh the pixel and consequently adjust the simplex so that its vertices better represent the homogeneous regions in natural images. The primary assumption of these methods is that the neighboring pixels of pure pixels are also pure or have some purity level that can boost the endmember extraction.

2) **No Pure Pixel Assumption:** If there are no pure pixels, virtual endmembers can be successfully located at the vertices if enough data points are available on the facets of the data simplex. Here, the main idea is to minimize the simplex volume subject to a convex combination of the endmembers containing the data points. For instance, the minimum volume enclosing simplex (MVES) [67], [68] seeks the minimum simplex, which encloses the data points. Maximum Volume Inscribed Ellipsoid (MVIE) [69] seeks the maximum volume ellipsoid contained within the convex hull of the data points. This volume is maximal when touching the facets of the data simplex. The contact points aid in estimating the virtual endmembers.

Alternatively, the minimum volume simplex analysis (MVSA) [70] and the simplex identification via variable splitting and augmented Lagrangian (SISAL) [24], [70] form a non-convex minimization (or a non-concave maximization)

problem which estimates endmembers by minimizing volume data simplex, subject to the ASC and ANC [24], [70]. Assuming $\mathbf{Q} = \mathbf{E}_r^{-1}$, SISAL solves

$$ \begin{gathered} \hat{\mathbf{Q}} = \arg\min_{\mathbf{Q}} -\log|\mathrm{det}(Q)| + \lambda \mathbf{1}_r^T \mathrm{hinge}(\mathbf{Q}\mathbf{Y}_r) \mathbf{1}_n \\ \text{s.t.} \quad \mathbf{1}_r^T \mathbf{Q} = \mathbf{1}_n^T \mathbf{Y}_r^{-1}. \end{gathered} \quad (25) $$

where the hinge function ($\mathrm{hinge}(x) = \max\{-x, 0\}$) is the regularization allowing violations to the non-negativity constraint and makes SISAL robust to noise. SISAL and MVSA solve the same problem however, SISAL is more efficient since it uses an alternating direction method of multipliers (ADMM) [71].

### C. Abundance Estimation

1) **Least Squares-based Approaches:** In supervised unmixing, when the endmember $\mathbf{E}$ are known, the abundances should be estimated. This step is also known as inversion, referring to the inverse problem of estimating abundances when the endmembers are known. The presence of noise, inevitable errors in endmember estimation/extraction, and the physical constraints make such an inverse problem very challenging.

Unconstrained least squares unmixing (UCLSU) via the orthogonal subspace projection was proposed for abundance estimation [72]. Non-negative constrained least squares unmixing (NCLSU) [73], [74] was proposed to estimate the abundances subjected to ANC. A weighted least square is suggested in [75] for estimating the abundances of multispectral remotely sensed data. There are several attempts to solve the least squares problem subjected to both ANC and ASC [75], [76]. In [75], quadratic programming was suggested for the constrained least squares. The first efficient algorithm was proposed in [22] and called fully constrained least squares unmixing (FCLSU). Later, fully constrained least squares unmixing by simplex projection was proposed [77] to derive a more efficient algorithm. A recursive algorithm minimizes the least square by performing orthogonal projections while holding the ASC and ANC. However, we should note that with the advances in graphical processing units (GPU), FCLSU can be efficiently solved using convex optimization techniques. Therefore, FCLSU is the most widely used method for abundance estimation given by:

$$ \begin{gathered} \hat{\mathbf{A}} = \arg\min_{\mathbf{A}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 \\ \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T. \end{gathered} \quad (26) $$

where the $||.||_F$ denotes the Frobenius norm. Problem (26) is convex and can be solved using any convex optimization solver. To promote spatial regularization, one can solve a penalized least square for abundance estimation given by

$$ \begin{gathered} \hat{\mathbf{A}} = \arg\min_{\mathbf{A}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \lambda\phi(\mathbf{A}) \\ \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \end{gathered} \quad (27) $$

where the first term is the fidelity term. The penalty term is defined based on function $\phi$ applied to the unknown abundances. $\phi$ can be a total variation or sparsity-based function
simplex inside the data. Simplex volume maximization (SiVM) [23] extracts the endmembers by iteratively maximizing the simplex volume using

$$ \arg \max_{\mathbf{E}_r} V(\mathbf{E}_r) = \arg \max_{\mathbf{E}_r} \sqrt{\frac{(-\mathbf{1})^r \cdot \mathrm{cmd}(\mathbf{E}_r)}{2^{r-1}(\mathbf{r}-1)!}}, \quad (24) $$

where cmd is the Cayley–Menger determinant:

$$ \mathrm{cmd}(\mathbf{E}_r) = \det \begin{bmatrix} 0 & 1 & 1 & 1 & \cdots & 1 \\ 1 & 0 & d_{1,2}^2 & d_{1,3}^2 & \cdots & d_{1,r}^2 \\ 1 & d_{2,1}^2 & 0 & d_{2,3}^2 & \cdots & d_{2,r}^2 \\ 1 & d_{3,1}^2 & d_{3,2}^2 & 0 & \cdots & d_{3,r}^2 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & d_{r,1}^2 & d_{r,2}^2 & d_{r,3}^2 & \cdots & 0 \end{bmatrix}, $$

and $d_{i,j}$ is the Euclidean distance between endmembers $e_i$ and $e_j$.

• **Multiple Pure Pixels and Endmember Bundles:** The multiple endmembers were also referred to as the endmember bundles. The endmember bundle and bundle unmixing were proposed in [58]–[60]. PCA was used to identify the (projected) endmembers using the extremes. More than one extreme was selected for materials. The bundle unmixing was formed using linear programming to determine the minimum, average, and maximum of the fractional abundances induced by endmember variability. In [61], endmember bundles were extracted from the observed data. Endmember variability is modeled using the endmember bundle, which will be discussed in more detail in semisupervised unmixing.

Some techniques incorporate spatial information for extracting the endmembers [62]–[65]. These methods often utilize information from neighboring pixels in the endmember extraction. For instance, in [66], the similarity of the pixel and its neighbors was used to weigh the pixel and consequently adjust the simplex so that its vertices better represent the homogeneous regions in natural images. The primary assumption of these methods is that the neighboring pixels of pure pixels are also pure or have some purity level that can boost the endmember extraction.

2) **No Pure Pixel Assumption:** If there are no pure pixels, virtual endmembers can be successfully located at the vertices if enough data points are available on the facets of the data simplex. Here, the main idea is to minimize the simplex volume subject to a convex combination of the endmembers containing the data points. For instance, the minimum volume enclosing simplex (MVES) [67], [68] seeks the minimum simplex, which encloses the data points. Maximum Volume Inscribed Ellipsoid (MVIE) [69] seeks the maximum volume ellipsoid contained within the convex hull of the data points. This volume is maximal when touching the facets of the data simplex. The contact points aid in estimating the virtual endmembers.

Alternatively, the minimum volume simplex analysis (MVSA) [70] and the simplex identification via variable splitting and augmented Lagrangian (SISAL) [24], [70] form a non-convex minimization (or a non-concave maximization)

problem which estimates endmembers by minimizing volume data simplex, subject to the ASC and ANC [24], [70]. Assuming $\mathbf{Q} = \mathbf{E}_r^{-1}$, SISAL solves

$$ \begin{gathered} \hat{\mathbf{Q}} = \arg \min_{\mathbf{Q}} -\log|\det(Q)| + \lambda \mathbf{1}_r^T \mathrm{hinge}(\mathbf{Q} \mathbf{Y}_r) \mathbf{1}_n \\ \text{s.t.} \quad \mathbf{1}_r^T \mathbf{Q} = \mathbf{1}_n^T \mathbf{Y}_r^{-1}. \end{gathered} \quad (25) $$

where the hinge function ($\mathrm{hinge}(x) = \max\{-x, 0\}$) is the regularization allowing violations to the non-negativity constraint and makes SISAL robust to noise. SISAL and MVSA solve the same problem however, SISAL is more efficient since it uses an alternating direction method of multipliers (ADMM) [71].

### C. Abundance Estimation

1) **Least Squares-based Approaches:** In supervised unmixing, when the endmember $\mathbf{E}$ are known, the abundances should be estimated. This step is also known as inversion, referring to the inverse problem of estimating abundances when the endmembers are known. The presence of noise, inevitable errors in endmember estimation/extraction, and the physical constraints make such an inverse problem very challenging.

Unconstrained least squares unmixing (UCLSU) via the orthogonal subspace projection was proposed for abundance estimation [72]. Non-negative constrained least squares unmixing (NCLSU) [73], [74] was proposed to estimate the abundances subjected to ANC. A weighted least square is suggested in [75] for estimating the abundances of multispectral remotely sensed data. There are several attempts to solve the least squares problem subjected to both ANC and ASC [75], [76]. In [75], quadratic programming was suggested for the constrained least squares. The first efficient algorithm was proposed in [22] and called fully constrained least squares unmixing (FCLSU). Later, fully constrained least squares unmixing by simplex projection was proposed [77] to derive a more efficient algorithm. A recursive algorithm minimizes the least square by performing orthogonal projections while holding the ASC and ANC. However, we should note that with the advances in graphical processing units (GPU), FCLSU can be efficiently solved using convex optimization techniques. Therefore, FCLSU is the most widely used method for abundance estimation given by:

$$ \begin{gathered} \hat{\mathbf{A}} = \arg \min_{\mathbf{A}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 \\ \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T. \end{gathered} \quad (26) $$

where the $||.||_F$ denotes the Frobenius norm. Problem (26) is convex and can be solved using any convex optimization solver. To promote spatial regularization, one can solve a penalized least square for abundance estimation given by

$$ \begin{gathered} \hat{\mathbf{A}} = \arg \min_{\mathbf{A}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \lambda \phi(\mathbf{A}) \\ \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \end{gathered} \quad (27) $$

where the first term is the fidelity term. The penalty term is defined based on function $\phi$ applied to the unknown abundances. $\phi$ can be a total variation or sparsity-based function
that usually incorporates spatial information. $\lambda$ controls the trade-off between the fidelity and penalty terms. Problem (27) is convex ($\mathbf{E}$ is known) as long as the penalty term is convex and, therefore, it can be solved using any convex optimization, least squares, or quadratic programming solver. Problem (27) will lead to better estimation in terms of RMSE compared to (26) for some values of $\lambda$. However, the selection of optimum $\lambda$ is not trivial.

2) **Shallow/Deep Neural Networks:** In supervised unmixing networks, the endmembers or the decoder weights are fixed. Therefore, a regularized loss function can be given by

$$ \mathcal{L}_{reg}(\mathbf{y}, \hat{\mathbf{y}}) = \mathcal{L}(\mathbf{y}, \hat{\mathbf{y}}) + \lambda J(\mathbf{a}, \mathbf{W}_{en}). \quad (28) $$

The regularizer $J$ can be a combination of a Frobenius norm on $\mathbf{W}_{en}$ with sparsity norms such as $\ell_1$ norm on $\mathbf{a}$. Autoencoders only utilize spectral information. Therefore, convolutional autoencoders were suggested for spectral unmixing to use both spatial and spectral information. Deep convolutional autoencoders were used for supervised hyperspectral unmixing in [78], exploiting 3D convolutional filters.

In [25], unmixing using deep image prior was proposed (UnDIP). Selection of the prior $\phi$ can be data dependent. Inspired by deep image prior [79], [80], UnDIP shows that the selection of the image prior $\phi$ in the minimization problem (27) can be shifted into optimizing a deep network's parameters. Then we have,

$$ \hat{\theta} = \arg \min_{\theta} \frac{1}{2} \|\mathbf{Y} - \mathbf{E} f_{\theta}(\mathbf{Z})\|_{F}^{2} \quad \text{s.t.} \quad \hat{\mathbf{A}} = f_{\theta}(\mathbf{Z}). \quad (29) $$

$\mathbf{Z}$ is the network input fixed throughout the learning. $f_{\theta}$ is a deep network with parameters $\theta$ representing the weights and biases. UnDIP uses SiVM to extract $\mathbf{E}$. Therefore, it is suitable for pure pixel-based scenarios. UnDIP uses an encoder-decoder network for $f_{\theta}$ and (29) as the loss function and a softmax layer to enforce ASC and ANC. In [81], a two-stream network was proposed. The network contains an encoder in parallel with an autoencoder. First, bundle endmembers were extracted using VCA, and the corresponding abundances were estimated using the linear mixing model. Then, the encoder was trained based on the bundle endmembers and the abundances. The encoder shares the weights with the encoder of the autoencoder, whose decoder is nonlinear. A similar idea was also used in [82].

## V. UNSUPERVISED (BLIND) UNMIXING

In supervised unmixing, the processing chain is sequential. Usually, the estimation of abundance does not affect the estimation of endmembers due to the order in the processing chain. In unsupervised unmixing, we assume that both endmember and abundances are unknown and blind unmixing endmember and abundances are estimated simultaneously. We consider three major paradigms; 1) Least Squares-based Approaches 2) Shallow/Deep learning-based Approaches, and 3) Statistical-Based Approaches. Due to the inherent nonconvexity of blind unmixing methods, they are often vulnerable to initialization, and therefore they are always initialized using a geometrical endmember extraction approach.

## A. Least Squares-based Approaches

Least squares-based approaches usually solve a minimization problem with respect to both $\mathbf{E}$ and $\mathbf{A}$ in the form of

$$ \begin{gathered} (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \lambda \phi(\mathbf{A}) + \beta \psi(\mathbf{E}) \\ \text{s.t.} \mathbf{A} \ge 0, \mathbf{1}_{r}^{T} \mathbf{A} = \mathbf{1}_{n}^{T}, 0 \le \mathbf{E} \le 1 \end{gathered} \quad (30) $$

The main difference between the unmixing approaches in this group is the selection of the penalty functions $\phi$ and $\psi$. $\phi$ is often selected for incorporating spatial information and capturing the spatial correlation. The common choices for $\phi$ are TV and sparsity-promoting norms. $\psi$ is usually a geometrical regularizer for incorporating geometrical information.

1) **Selection of a Regularizer:** The common choices for $\psi$ are simplex volume terms and norms with simplex volume flavor, such as TV [83]. Total variation penalty enforces the data simplex to have minimum volume [83]

$$ \begin{gathered} (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \beta \mathrm{TV}_{e}(\mathbf{E}) \\ \text{s.t.} \mathbf{A} \ge 0, \mathbf{1}_{r}^{T} \mathbf{A} = \mathbf{1}_{n}^{T} \end{gathered} \quad (31) $$

where $\mathrm{TV}_e(\mathbf{E})$ is a TV function given by

$$ \mathrm{TV}_e(\mathbf{E}) = \sum_{i,j=1}^{r} \|\mathbf{e}_i - \mathbf{e}_j\|_2^2 = \|\mathbf{E}(\mathbf{I}_r - \frac{1}{r}\mathbf{1}_r\mathbf{1}_r^T)\|_F^2. \quad (32) $$

See Appendix A for more details. We should mention that the method proposed in [83] solves the problem after subspace projection using MNF.

Minimum volume-constrained nonnegative matrix factorization (MVC-NMF) was proposed in [84]. MVC-NMF suggests a minimum volume penalty for the constraint nonnegative least squares given by

$$ \begin{gathered} (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \beta \mathrm{Vol}^{2}(\mathbf{E}_{r}) \\ \text{s.t.} \mathbf{A} \ge 0, \mathbf{1}_{r}^{T} \mathbf{A} = \mathbf{1}_{n}^{T}, \mathbf{E} \ge 0. \end{gathered} \quad (33) $$

Collaborative nonnegative matrix factorization (CoMNF) [85] suggests a sparsity-promoting term on abundances and MV on the simplex given by

$$ \begin{gathered} (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \beta \|\mathbf{E} - \mathbf{m}\mathbf{1}_{r}^{T}\|_{F}^{2} \\ + \lambda \sum_{i=1}^{r} \|\mathbf{a}_{(i)}\|_{2}^{q} \quad \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_{r}^{T} \mathbf{A} = \mathbf{1}_{n}^{T}, \end{gathered} \quad (34) $$

where, $0 \le q \le 1$, $\mathbf{m}$ contains the mean values of the spectral pixels, i.e., $\mathbf{m} = \frac{1}{n}\mathbf{Y}\mathbf{1}_{n}$, and $\mathbf{a}_{(i)}$ is the $i$th row of the matrix $\mathbf{A}$. This term pulls the endmembers toward the center of mass. CoNMF uses both spatial and MV (geometrical) regularizers and solves the problem by projecting the data into a subspace. In [86], Robust CoMNF (RCoNMF) was proposed, which utilizes a geometrical penalty that minimizes the distances between the endmembers to be estimated and the boundary pixels ($\mathbf{P} \in \mathbb{R}^{p \times r}$). The main assumption is that
that usually incorporates spatial information. $\lambda$ controls the trade-off between the fidelity and penalty terms. Problem (27) is convex ($\mathbf{E}$ is known) as long as the penalty term is convex and, therefore, it can be solved using any convex optimization, least squares, or quadratic programming solver. Problem (27) will lead to better estimation in terms of RMSE compared to (26) for some values of $\lambda$. However, the selection of optimum $\lambda$ is not trivial.

2) **Shallow/Deep Neural Networks:** In supervised unmixing networks, the endmembers or the decoder weights are fixed. Therefore, a regularized loss function can be given by

$$ \mathcal{L}_{reg}(\mathbf{y}, \hat{\mathbf{y}}) = \mathcal{L}(\mathbf{y}, \hat{\mathbf{y}}) + \lambda J(\mathbf{a}, \mathbf{W}_{en}). \quad (28) $$

The regularizer $J$ can be a combination of a Frobenius norm on $\mathbf{W}_{en}$ with sparsity norms such as $\ell_1$ norm on $\mathbf{a}$. Autoencoders only utilize spectral information. Therefore, convolutional autoencoders were suggested for spectral unmixing to use both spatial and spectral information. Deep convolutional autoencoders were used for supervised hyperspectral unmixing in [78], exploiting 3D convolutional filters.

In [25], unmixing using deep image prior was proposed (UnDIP). Selection of the prior $\phi$ can be data dependent. Inspired by deep image prior [79], [80], UnDIP shows that the selection of the image prior $\phi$ in the minimization problem (27) can be shifted into optimizing a deep network's parameters. Then we have,

$$ \hat{\theta} = \arg \min_{\theta} \frac{1}{2} \|\mathbf{Y} - \mathbf{E} f_{\theta}(\mathbf{Z})\|_{F}^{2} \quad \text{s.t.} \quad \hat{\mathbf{A}} = f_{\theta}(\mathbf{Z}). \quad (29) $$

$\mathbf{Z}$ is the network input fixed throughout the learning. $f_{\theta}$ is a deep network with parameters $\theta$ representing the weights and biases. UnDIP uses SiVM to extract $\mathbf{E}$. Therefore, it is suitable for pure pixel-based scenarios. UnDIP uses an encoder-decoder network for $f_{\theta}$ and (29) as the loss function and a softmax layer to enforce ASC and ANC. In [81], a two-stream network was proposed. The network contains an encoder in parallel with an autoencoder. First, bundle endmembers were extracted using VCA, and the corresponding abundances were estimated using the linear mixing model. Then, the encoder was trained based on the bundle endmembers and the abundances. The encoder shares the weights with the encoder of the autoencoder, whose decoder is nonlinear. A similar idea was also used in [82].

## V. UNSUPERVISED (BLIND) UNMIXING

In supervised unmixing, the processing chain is sequential. Usually, the estimation of abundance does not affect the estimation of endmembers due to the order in the processing chain. In unsupervised unmixing, we assume that both endmember and abundances are unknown and blind unmixing endmember and abundances are estimated simultaneously. We consider three major paradigms; 1) Least Squares-based Approaches 2) Shallow/Deep learning-based Approaches, and 3) Statistical-Based Approaches. Due to the inherent nonconvexity of blind unmixing methods, they are often vulnerable to initialization, and therefore they are always initialized using a geometrical endmember extraction approach.

## A. Least Squares-based Approaches

Least squares-based approaches usually solve a minimization problem with respect to both $\mathbf{E}$ and $\mathbf{A}$ in the form of

$$ \begin{gathered} (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \lambda \phi(\mathbf{A}) + \beta \psi(\mathbf{E}) \\ \text{s.t.} \mathbf{A} \ge 0, \mathbf{1}_{r}^{T} \mathbf{A} = \mathbf{1}_{n}^{T}, 0 \le \mathbf{E} \le 1 \end{gathered} \quad (30) $$

The main difference between the unmixing approaches in this group is the selection of the penalty functions $\phi$ and $\psi$. $\phi$ is often selected for incorporating spatial information and capturing the spatial correlation. The common choices for $\phi$ are TV and sparsity-promoting norms. $\psi$ is usually a geometrical regularizer for incorporating geometrical information.

1) **Selection of a Regularizer:** The common choices for $\psi$ are simplex volume terms and norms with simplex volume flavor, such as TV [83]. Total variation penalty enforces the data simplex to have minimum volume [83]

$$ \begin{gathered} (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \beta \mathrm{TV}_{e}(\mathbf{E}) \\ \text{s.t.} \mathbf{A} \ge 0, \mathbf{1}_{r}^{T} \mathbf{A} = \mathbf{1}_{n}^{T} \end{gathered} \quad (31) $$

where $\mathrm{TV}_e(\mathbf{E})$ is a TV function given by

$$ \mathrm{TV}_e(\mathbf{E}) = \sum_{i,j=1}^{r} \|\mathbf{e}_i - \mathbf{e}_j\|_2^2 = \|\mathbf{E}(\mathbf{I}_r - \frac{1}{r}\mathbf{1}_r\mathbf{1}_r^T)\|_F^2. \quad (32) $$

See Appendix A for more details. We should mention that the method proposed in [83] solves the problem after subspace projection using MNF.

Minimum volume-constrained nonnegative matrix factorization (MVC-NMF) was proposed in [84]. MVC-NMF suggests a minimum volume penalty for the constraint nonnegative least squares given by

$$ \begin{gathered} (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \beta \mathrm{Vol}^{2}(\mathbf{E}_{r}) \\ \text{s.t.} \mathbf{A} \ge 0, \mathbf{1}_{r}^{T} \mathbf{A} = \mathbf{1}_{n}^{T}, \mathbf{E} \ge 0. \end{gathered} \quad (33) $$

Collaborative nonnegative matrix factorization (CoMNF) [85] suggests a sparsity-promoting term on abundances and MV on the simplex given by

$$ \begin{gathered} (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \beta \|\mathbf{E} - \mathbf{m}\mathbf{1}_{r}^{T}\|_{F}^{2} \\ + \lambda \sum_{i=1}^{r} \|\mathbf{a}_{(i)}\|_{2}^{q} \quad \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_{r}^{T} \mathbf{A} = \mathbf{1}_{n}^{T}, \end{gathered} \quad (34) $$

where, $0 \le q \le 1$, $\mathbf{m}$ contains the mean values of the spectral pixels, i.e., $\mathbf{m} = \frac{1}{n}\mathbf{Y}\mathbf{1}_{n}$, and $\mathbf{a}_{(i)}$ is the $i$th row of the matrix $\mathbf{A}$. This term pulls the endmembers toward the center of mass. CoNMF uses both spatial and MV (geometrical) regularizers and solves the problem by projecting the data into a subspace. In [86], Robust CoMNF (RCoNMF) was proposed, which utilizes a geometrical penalty that minimizes the distances between the endmembers to be estimated and the boundary pixels ($\mathbf{P} \in \mathbb{R}^{p \times r}$). The main assumption is that
the endmembers are close to the extremes of the data simplex (so-called boundary pixels). Hence, RCoNMF solves

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \beta \|\mathbf{E} - \mathbf{P}\|_F^2 \\ + \lambda \sum_{i=1}^r \|\mathbf{a}_{(i)}\|_2 \quad \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \qquad (35) $$

where the boundary pixels, **P**, can be obtained using pure pixel-based endmember extraction. In [86], VCA was used to estimate **P**. In [87], $\ell_{0.5}$-NMF was proposed which enforces $\ell_q = \sum_{i=1}^r \sum_{j=1}^n |a_{ij}|^{0.5}$ sparsity penalty ($q = 1/2$) on the abundances. A multiplicative update rule [88], [89] was proposed to solve the $\ell_{0.5}$-NMF problem.

In such minimization problems, selecting the regulation parameters is not trivial. In [38], Nonnegative Matrix Factorization-Quadratic Minimum Volume (NMF-QMV) was proposed. NMF-QMV unifies the three geometrical regularizers mentioned above (i.e., TV, center, and boundary) in a single framework. NMF-QMV proposes the unified cost function using the MV-regularizers given by

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \beta \|\mathbf{E}\mathbf{G} - \mathbf{O}\|_F^2 \\ \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \qquad (36) $$

where the selection of **G** and **O** are given in Table II. We emphasize that those terms do not compute the volume of a simplex. They are convex penalties with volume flavor. A parameter selection technique was proposed in 36 by searching for an optimum $\beta$ in a large interval. $\beta$ optimum is when the estimated endmember simplex boundary is close to the data points boundary. Note that the spatial regularizer was ignored in (36), simplifying the selection technique. We

<table><caption>TABLE II
SELECTION OF **G** AND **O** FOR THE UNIFIED GEOMETRICAL REGULARIZER.</caption><thead><tr><th></th><th>TV</th><th>Center</th><th>Boundry</th></tr></thead><tbody><tr><th><b>G</b></th><td><b>I</b><sub>r</sub> - <sup>1</sup>&frasl;<sub>r</sub><b>1</b><sub>r</sub><b>1</b><sub>r</sub><sup>T</sup></td><td><b>I</b><sub>r</sub></td><td><b>I</b><sub>r</sub></td></tr><tr><th><b>O</b></th><td><b>0</b></td><td><b>m</b><b>1</b><sub>r</sub></td><td>Extracted Extremes</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Fig. 12. The intuition of how the three geometrical regularizers enforce the endmembers using (a) TV, (b) center, and (c) boundary terms to change an initial simplex (blue) towards the endmember simplex (red). The black arrows show the direction of the force.</figcaption></figure>

should note that some blind unmixing techniques [38], [83], [85], [86], solve the problems in a subspace. This has two main advantages: 1) reducing the computations, and 2) reducing noise and removing outliers. Hence, the endmembers are assumed to belong to an $r-1$ dimensional affine set, which best represents the data using an orthogonal projection. Note

that the abundances do not change and can be estimated in a subspace. However, in real-world scenarios, the estimated endmembers in a subspace after reprojection into the original space may have negative values due to noise and nonlinearity since they cannot be bounded between zero and one in a subspace.

In [90], the spectral roughness penalty was proposed for imposing the smoothness to the endmembers, and the $\ell_q$ norm was applied to the abundances to enforce sparsity on the abundances. The proposed cost function is given by

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \frac{\beta}{2} \|\mathbf{R}\mathbf{E}\|_F^2 + \lambda \sum_{i=1}^n \|\mathbf{a}_i\|_q \\ + \frac{\alpha}{2} \|\mathbf{1}_r^T \mathbf{A} - \mathbf{1}_n^T\|_2^2 \quad \text{s.t.} \quad \mathbf{A} \ge 0, \qquad (37) $$

where $0 \le q \le 1$ and $\mathbf{R} \in \mathbb{R}^{(p-1) \times p}$ is the first-order difference matrix given by

$$ \mathbf{R} = \begin{bmatrix} -1 & 1 & 0 & \cdots & 0 & 0 \\ 0 & -1 & 1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & -1 & 1 \end{bmatrix}. \qquad (38) $$

We should note that the roughness penalty does not incorporate geometrical information, and ASC constraint was relaxed and added to the cost function as a penalty. This makes the algorithm robust to the noise however, it adds one extra regularization parameter, $\alpha$, which together with $\beta$ and $\lambda$ become very challenging to select, particularly in real-world applications. In [91], the (2 dimensional) isotropic total variation was applied to abundances to promote piece-wise smoothness spatially. The cost function is given by

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \lambda_1 \sum_{i=1}^r TV_i(\mathbf{a}_{(i)}) + \\ \lambda_2 \sum_{i=1}^r \|\mathbf{a}_{(i)}\|_q \quad \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{E} \ge 0, \|\mathbf{e}_{(i)}\|_2^2 = 1, i = 1, 2, \dots, r, \qquad (39) $$

where $0 \le q \le 1$ and $TV_i$ is the nonisotropic TV function given by

$$ TV_i(\mathbf{x}) = \left\| \sqrt{(\mathbf{H}_h(\mathbf{x}))^2 + (\mathbf{H}_v(\mathbf{x}))^2} \right\|_1. \qquad (40) $$

$\mathbf{H}_v = \mathbf{I} \otimes \mathbf{R}$ and $\mathbf{H}_h = \mathbf{R} \otimes \mathbf{I}$ are the matrix operators to calculate the first order vertical and horizontal differences, respectively, for the vectorized image $\mathbf{x}$. We should note that ASC was ignored in (39), and the endmember norm constraint (ENC) was applied to enforce a unit norm on the endmembers. ENC captures the spectral variability but has a big disadvantage, which does not let the algorithm distinguish the materials with similar endmembers with some scaling factors, such as tree and grass. In [92], NMF subject to ASC was suggested using a reweighted $\ell_1$ norm and the total variation penalty applied to the abundances. Many NMF-based blind unmixing approaches were proposed using the form of minimization (30) [93]–[96].
the endmembers are close to the extremes of the data simplex (so-called boundary pixels). Hence, RCoNMF solves

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \beta \|\mathbf{E} - \mathbf{P}\|_F^2 \\ + \lambda \sum_{i=1}^r \|\mathbf{a}_{(i)}\|_2 \quad \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \qquad (35) $$

where the boundary pixels, **P**, can be obtained using pure pixel-based endmember extraction. In [86], VCA was used to estimate **P**. In [87], $\ell_{0.5}$-NMF was proposed which enforces $\ell_q = \sum_{i=1}^r \sum_{j=1}^n |a_{ij}|^{0.5}$ sparsity penalty ($q = 1/2$) on the abundances. A multiplicative update rule [88], [89] was proposed to solve the $\ell_{0.5}$-NMF problem.

In such minimization problems, selecting the regulation parameters is not trivial. In [38], Nonnegative Matrix Factorization-Quadratic Minimum Volume (NMF-QMV) was proposed. NMF-QMV unifies the three geometrical regularizers mentioned above (i.e., TV, center, and boundary) in a single framework. NMF-QMV proposes the unified cost function using the MV-regularizers given by

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \beta \|\mathbf{E}\mathbf{G} - \mathbf{O}\|_F^2 \\ \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \qquad (36) $$

where the selection of **G** and **O** are given in Table II. We emphasize that those terms do not compute the volume of a simplex. They are convex penalties with volume flavor. A parameter selection technique was proposed in 36 by searching for an optimum $\beta$ in a large interval. $\beta$ optimum is when the estimated endmember simplex boundary is close to the data points boundary. Note that the spatial regularizer was ignored in (36), simplifying the selection technique. We

<table><caption>TABLE II
SELECTION OF **G** AND **O** FOR THE UNIFIED GEOMETRICAL REGULARIZER.</caption><thead><tr><th></th><th>TV</th><th>Center</th><th>Boundry</th></tr></thead><tbody><tr><th><b>G</b></th><td><b>I</b><sub>r</sub> - <sup>1</sup>&frasl;<sub>r</sub><b>1</b><sub>r</sub><b>1</b><sub>r</sub><sup>T</sup></td><td><b>I</b><sub>r</sub></td><td><b>I</b><sub>r</sub></td></tr><tr><th><b>O</b></th><td><b>0</b></td><td><b>m</b><b>1</b><sub>r</sub></td><td>Extracted Extremes</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Fig. 12. The intuition of how the three geometrical regularizers enforce the endmembers using (a) TV, (b) center, and (c) boundary terms to change an initial simplex (blue) towards the endmember simplex (red). The black arrows show the direction of the force.</figcaption></figure>

should note that some blind unmixing techniques [38], [83], [85], [86], solve the problems in a subspace. This has two main advantages: 1) reducing the computations, and 2) reducing noise and removing outliers. Hence, the endmembers are assumed to belong to an $r-1$ dimensional affine set, which best represents the data using an orthogonal projection. Note

that the abundances do not change and can be estimated in a subspace. However, in real-world scenarios, the estimated endmembers in a subspace after reprojection into the original space may have negative values due to noise and nonlinearity since they cannot be bounded between zero and one in a subspace.

In [90], the spectral roughness penalty was proposed for imposing the smoothness to the endmembers, and the $\ell_q$ norm was applied to the abundances to enforce sparsity on the abundances. The proposed cost function is given by

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \frac{\beta}{2} \|\mathbf{R}\mathbf{E}\|_F^2 + \lambda \sum_{i=1}^n \|\mathbf{a}_i\|_q \\ + \frac{\alpha}{2} \|\mathbf{1}_r^T \mathbf{A} - \mathbf{1}_n^T\|_2^2 \quad \text{s.t.} \quad \mathbf{A} \ge 0, \qquad (37) $$

where $0 \le q \le 1$ and $\mathbf{R} \in \mathbb{R}^{(p-1) \times p}$ is the first-order difference matrix given by

$$ \mathbf{R} = \begin{bmatrix} -1 & 1 & 0 & \cdots & 0 & 0 \\ 0 & -1 & 1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & -1 & 1 \end{bmatrix}. \qquad (38) $$

We should note that the roughness penalty does not incorporate geometrical information, and ASC constraint was relaxed and added to the cost function as a penalty. This makes the algorithm robust to the noise however, it adds one extra regularization parameter, $\alpha$, which together with $\beta$ and $\lambda$ become very challenging to select, particularly in real-world applications. In [91], the (2 dimensional) isotropic total variation was applied to abundances to promote piece-wise smoothness spatially. The cost function is given by

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \lambda_1 \sum_{i=1}^r TV_i(\mathbf{a}_{(i)}) + \\ \lambda_2 \sum_{i=1}^r \|\mathbf{a}_{(i)}\|_q \quad \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{E} \ge 0, \|\mathbf{e}_{(i)}\|_2^2 = 1, i = 1, 2, \dots, r, \qquad (39) $$

where $0 \le q \le 1$ and $TV_i$ is the nonisotropic TV function given by

$$ TV_i(\mathbf{x}) = \left\| \sqrt{(\mathbf{H}_h(\mathbf{x}))^2 + (\mathbf{H}_v(\mathbf{x}))^2} \right\|_1. \qquad (40) $$

$\mathbf{H}_v = \mathbf{I} \otimes \mathbf{R}$ and $\mathbf{H}_h = \mathbf{R} \otimes \mathbf{I}$ are the matrix operators to calculate the first order vertical and horizontal differences, respectively, for the vectorized image $\mathbf{x}$. We should note that ASC was ignored in (39), and the endmember norm constraint (ENC) was applied to enforce a unit norm on the endmembers. ENC captures the spectral variability but has a big disadvantage, which does not let the algorithm distinguish the materials with similar endmembers with some scaling factors, such as tree and grass. In [92], NMF subject to ASC was suggested using a reweighted $\ell_1$ norm and the total variation penalty applied to the abundances. Many NMF-based blind unmixing approaches were proposed using the form of minimization (30) [93]–[96].
A plug-and-play (PnP) priors framework was proposed in [97] for blind unmixing. Problem (30), omitting the penalty and constraint on $E$, was solved using an ADMM. The penalized least squares step in the ADMM can be solved using a denoiser. Several denoisers, including a deep network, were suggested to solve that step. Additionally, the penalty term was considered not only applied to the abundances but also to the reconstructed data, i.e., $\phi(EA)$.

2) **Archetypal Analysis Models:** Several archetypal analysis (AA)-based [98] algorithms were proposed for blind unmixing. Let's assume that the endmembers are convex combinations of the spectral pixels, then using AA we have

$$ (\hat{\mathbf{A}}, \hat{\mathbf{B}}) = \arg \min_{\mathbf{A}, \mathbf{B}} \frac{1}{2} \|\mathbf{Y} - \mathbf{Y}\mathbf{A}\|_F^2 \\ \text{s.t.} \quad \mathbf{B} \ge 0, \mathbf{1}_n^T \mathbf{B} = \mathbf{1}_r^T, \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \quad (41) $$

where $\mathbf{B} \in \mathbb{R}^{n \times r}$ and the columns of $\mathbf{B}$ belong to the simplex $\Delta_n$. In [35], Entropic Descent Archetypal Analysis (EDAA) was proposed. EDAA uses an entropic descent algorithm to solve the AA problem (41). A model selection technique suggested where the coherence $\mu = \max_{i \neq j} e_i^T e_j$ and $\ell_1$ residual $\|\mathbf{Y} - \hat{\mathbf{E}}\hat{\mathbf{A}}\|_1$ are minimized. This makes the solution robust to the initialization of $\mathbf{B}$. The convex combination constraint on the endmembers is very strict. Therefore, in [99], a variation of AA suggested which uses a relaxation factor $\delta$ and an $r$ dimensional scaling vector $\mathbf{s} \in \mathbb{R}^r$

$$ (\hat{\mathbf{A}}, \hat{\mathbf{B}}, \hat{\mathbf{S}}) = \arg \min_{\mathbf{A}, \mathbf{B}, \mathbf{S}} \frac{1}{2} \|\mathbf{Y} - \mathbf{Y}\mathbf{B}(\mathbf{S} \odot \mathbf{I})\mathbf{A}\|_F^2 \text{ s.t.} \\ \mathbf{B} \ge 0, \mathbf{1}_n^T \mathbf{B} = \mathbf{1}_r^T, \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, 1 - \delta \le \mathbf{s} \le 1 + \delta \quad (42) $$

where $\mathbf{S} = \mathbf{s}\mathbf{1}^T$ and $\odot$ denotes the Hadamard product. Additionally, Kernel AA (KAA) was suggested to take into account the endmember variability by multiple endmember extraction. The proposed KAA was computationally expensive, and therefore, in [100], the Nyström method is used to construct a low-rank approximation of the high-dimensional kernel matrix using K-means.

An $\ell_1$ penalized variation of AA was proposed in [101] given by

$$ (\hat{\mathbf{A}}, \hat{\mathbf{B}}, \hat{\mathbf{S}}) = \arg \min_{\mathbf{A}, \mathbf{B}, \mathbf{S}} \frac{1}{2} \|\mathbf{Y} - \mathbf{P}\mathbf{B}(\mathbf{S} \odot \mathbf{I})\mathbf{A}\|_F^2 + \|\mathbf{A}\|_1 \quad \text{s.t.} \\ \mathbf{B} \ge 0, \mathbf{1}_n^T \mathbf{B} = \mathbf{1}_r^T, \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, 1 - \delta \le \mathbf{s} \le 1 + \delta \quad (43) $$

where the sparsity was enforced on the abundance and $\mathbf{P} \in \mathbb{R}^{n \times c}$, $c \gg r$ contains boundary pixels extracted by using PPI to cope with the endmember variability.

3) **Tensor Models:** Tensor decomposition (TD) was also suggested for blind unmixing motivated by the tensor structure of hyperspectral data, which can be represented as a three-dimensional array (tensor) with dimensions corresponding to spatial coordinates (row and columns) and spectral bands (channels) [102]. Hyperspectral data are often decomposed into nonnegative tensors. The interpretation of endmembers and abundances in a tensor decomposition makes the unmixing modeling challenging. The canonical polyadic decomposition (CPD) [102] and Tucker decomposition (or higher order SVD)

<figure><img src="image_14.png" /><figcaption>Fig. 13. The rank-$(L_i, L_i, 1)$ decomposition which is widely used for blind linear unmixing.</figcaption></figure>

[103] have been proposed for unmixing. However, they lack such an interpretation. Among tensor decompositions, matrix-vector third-order tensor decomposition (also known as rank-$(L_i, L_i, 1)$ block-term decomposition) is the most popular model for unmixing due to the interpretation of endmembers and abundances in the decomposition. It is a specific case of block term decompositions (BTDs) [104] which is generally defined as a sum of r rank-$(L_i, M_i, N_i)$ terms ($i = 1, 2, ..., r$). BTD can be assumed as a generalization of CPD and Tucker decomposition where it decomposes a tensor into a sum of component tensors (same as CPD), while each component tensor is factorized using a Tucker decomposition. Considering the hyperspectral data as a 3rd-order tensor $\mathcal{Y} \in \mathbb{R}^{n_1 \times n_2 \times p}$ (where $n_1 \times n_2 = n$), rank-$(L_i, L_i, 1)$ decomposition is written as

$$ \mathcal{Y} = \sum_{i=1}^{r} \mathbf{A}_i \odot \mathbf{e}_i = \sum_{i=1}^{r} \mathbf{K}_i \mathbf{B}_i^T \odot \mathbf{e}_i + \mathcal{N}, \mathbf{K}_i \ge 0, \mathbf{B}_i \ge 0, \mathbf{e}_i \ge 0 \quad (44) $$

where $\mathbf{A}_i$ denotes $i$th $n_1 \times n_2$ abundance (i.e., $vec^{-1}(\mathbf{a}_i)$) of rank $L_i$. Therefore, $\mathbf{A}_i$ can be written as a matrix product $\mathbf{A}_i = \mathbf{K}_i \mathbf{B}_i^T$ where $\mathbf{K}_i \in \mathbb{R}^{n_1 \times L_i}$ and $\mathbf{B}_i \in \mathbb{R}^{n_2 \times L_i}$ are full column rank matrices of rank $L_i$. The rank-$(L_i, L_i, 1)$ decomposition shown in Fig 13 is unique and a sufficient uniqueness condition is that the partitioned matrices $\mathbf{K} = [\mathbf{K}_1, \mathbf{K}_2, ..., \mathbf{K}_r]$ and $\mathbf{B} = [\mathbf{B}_1, \mathbf{B}_2, ..., \mathbf{B}_r]$ are full column rank and $E$ does not have any collinear columns [104]. It was shown that full column rank of $\mathbf{A}$ and $\mathbf{B}$ can be met if $\min(n_1, n_2) = \sum_{i=1}^r L_i$. Those conditions may be met in unmixing applications however they are not ensured. We should note that the abovementioned uniqueness conditions were discussed for BTD without nonnegativity constraints. Nevertheless, the final solution always depends on the initialization due to the nonconvexity in blind unmixing. Matrix-Vector Nonnegative Tensor Factorization (MVNTF) was proposed in [105] which solves

$$ (\hat{\mathbf{K}}, \hat{\mathbf{B}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{K}, \mathbf{B}, \mathbf{E}} \frac{1}{2} \|\mathcal{Y} - \sum_{i=1}^{r} \mathbf{K}_i \mathbf{B}_i^T \odot \mathbf{e}_i\|_F^2 + \\ \frac{\lambda}{2} \| \sum_{i=1}^{r} \mathbf{K}_i \mathbf{B}_i^T - \mathbf{1}_{n_1} \mathbf{1}_{n_2}^T \|_F^2 \quad \text{s.t.} \quad \mathbf{K} \ge 0, \mathbf{B} \ge 0, \mathbf{E} \ge 0 \quad (45) $$

where the Frobenius norm for a third-order tensor is given by $\mathcal{X} = \sqrt{\sum_i \sum_j \sum_k x_{ijk}^2}$ and the second term ensures the ASC. Several variations of (45) were proposed in the literature. In [106], sparse and low-rank NTF was proposed which enforces both $\ell_1$ sparsity and nuclear norm on abundances. This method was further improved in [107] by exploiting weighted norms
A plug-and-play (PnP) priors framework was proposed in [97] for blind unmixing. Problem (30), omitting the penalty and constraint on $E$, was solved using an ADMM. The penalized least squares step in the ADMM can be solved using a denoiser. Several denoisers, including a deep network, were suggested to solve that step. Additionally, the penalty term was considered not only applied to the abundances but also to the reconstructed data, i.e., $\phi(EA)$.

2) **Archetypal Analysis Models:** Several archetypal analysis (AA)-based [98] algorithms were proposed for blind unmixing. Let's assume that the endmembers are convex combinations of the spectral pixels, then using AA we have

$$ (\hat{\mathbf{A}}, \hat{\mathbf{B}}) = \arg \min_{\mathbf{A}, \mathbf{B}} \frac{1}{2} \|\mathbf{Y} - \mathbf{Y}\mathbf{A}\|_F^2 \\ \text{s.t.} \quad \mathbf{B} \ge 0, \mathbf{1}_n^T \mathbf{B} = \mathbf{1}_r^T, \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \quad (41) $$

where $\mathbf{B} \in \mathbb{R}^{n \times r}$ and the columns of $\mathbf{B}$ belong to the simplex $\Delta_n$. In [35], Entropic Descent Archetypal Analysis (EDAA) was proposed. EDAA uses an entropic descent algorithm to solve the AA problem (41). A model selection technique suggested where the coherence $\mu = \max_{i \neq j} e_i^T e_j$ and $\ell_1$ residual $\|\mathbf{Y} - \hat{\mathbf{E}}\hat{\mathbf{A}}\|_1$ are minimized. This makes the solution robust to the initialization of $\mathbf{B}$. The convex combination constraint on the endmembers is very strict. Therefore, in [99], a variation of AA suggested which uses a relaxation factor $\delta$ and an $r$ dimensional scaling vector $\mathbf{s} \in \mathbb{R}^r$

$$ (\hat{\mathbf{A}}, \hat{\mathbf{B}}, \hat{\mathbf{S}}) = \arg \min_{\mathbf{A}, \mathbf{B}, \mathbf{S}} \frac{1}{2} \|\mathbf{Y} - \mathbf{Y}\mathbf{B}(\mathbf{S} \odot \mathbf{I})\mathbf{A}\|_F^2 \text{ s.t.} \\ \mathbf{B} \ge 0, \mathbf{1}_n^T \mathbf{B} = \mathbf{1}_r^T, \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, 1 - \delta \le \mathbf{s} \le 1 + \delta \quad (42) $$

where $\mathbf{S} = \mathbf{s}\mathbf{1}^T$ and $\odot$ denotes the Hadamard product. Additionally, Kernel AA (KAA) was suggested to take into account the endmember variability by multiple endmember extraction. The proposed KAA was computationally expensive, and therefore, in [100], the Nyström method is used to construct a low-rank approximation of the high-dimensional kernel matrix using K-means.

An $\ell_1$ penalized variation of AA was proposed in [101] given by

$$ (\hat{\mathbf{A}}, \hat{\mathbf{B}}, \hat{\mathbf{S}}) = \arg \min_{\mathbf{A}, \mathbf{B}, \mathbf{S}} \frac{1}{2} \|\mathbf{Y} - \mathbf{P}\mathbf{B}(\mathbf{S} \odot \mathbf{I})\mathbf{A}\|_F^2 + \|\mathbf{A}\|_1 \quad \text{s.t.} \\ \mathbf{B} \ge 0, \mathbf{1}_n^T \mathbf{B} = \mathbf{1}_r^T, \mathbf{A} \ge 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, 1 - \delta \le \mathbf{s} \le 1 + \delta \quad (43) $$

where the sparsity was enforced on the abundance and $\mathbf{P} \in \mathbb{R}^{n \times c}$, $c \gg r$ contains boundary pixels extracted by using PPI to cope with the endmember variability.

3) **Tensor Models:** Tensor decomposition (TD) was also suggested for blind unmixing motivated by the tensor structure of hyperspectral data, which can be represented as a three-dimensional array (tensor) with dimensions corresponding to spatial coordinates (row and columns) and spectral bands (channels) [102]. Hyperspectral data are often decomposed into nonnegative tensors. The interpretation of endmembers and abundances in a tensor decomposition makes the unmixing modeling challenging. The canonical polyadic decomposition (CPD) [102] and Tucker decomposition (or higher order SVD)

<figure><img src="image_14.png" /><figcaption>Fig. 13. The rank-$(L_i, L_i, 1)$ decomposition which is widely used for blind linear unmixing.</figcaption></figure>

[103] have been proposed for unmixing. However, they lack such an interpretation. Among tensor decompositions, matrix-vector third-order tensor decomposition (also known as rank-$(L_i, L_i, 1)$ block-term decomposition) is the most popular model for unmixing due to the interpretation of endmembers and abundances in the decomposition. It is a specific case of block term decompositions (BTDs) [104] which is generally defined as a sum of r rank-$(L_i, M_i, N_i)$ terms ($i = 1, 2, ..., r$). BTD can be assumed as a generalization of CPD and Tucker decomposition where it decomposes a tensor into a sum of component tensors (same as CPD), while each component tensor is factorized using a Tucker decomposition. Considering the hyperspectral data as a 3rd-order tensor $\mathcal{Y} \in \mathbb{R}^{n_1 \times n_2 \times p}$ (where $n_1 \times n_2 = n$), rank-$(L_i, L_i, 1)$ decomposition is written as

$$ \mathcal{Y} = \sum_{i=1}^{r} \mathbf{A}_i \odot \mathbf{e}_i = \sum_{i=1}^{r} \mathbf{K}_i \mathbf{B}_i^T \odot \mathbf{e}_i + \mathcal{N}, \mathbf{K}_i \ge 0, \mathbf{B}_i \ge 0, \mathbf{e}_i \ge 0 \quad (44) $$

where $\mathbf{A}_i$ denotes $i$th $n_1 \times n_2$ abundance (i.e., $\text{vec}^{-1}(\mathbf{a}_i)$) of rank $L_i$. Therefore, $\mathbf{A}_i$ can be written as a matrix product $\mathbf{A}_i = \mathbf{K}_i \mathbf{B}_i^T$ where $\mathbf{K}_i \in \mathbb{R}^{n_1 \times L_i}$ and $\mathbf{B}_i \in \mathbb{R}^{n_2 \times L_i}$ are full column rank matrices of rank $L_i$. The rank-$(L_i, L_i, 1)$ decomposition shown in Fig 13 is unique and a sufficient uniqueness condition is that the partitioned matrices $\mathbf{K} = [\mathbf{K}_1, \mathbf{K}_2, ..., \mathbf{K}_r]$ and $\mathbf{B} = [\mathbf{B}_1, \mathbf{B}_2, ..., \mathbf{B}_r]$ are full column rank and $E$ does not have any collinear columns [104]. It was shown that full column rank of $\mathbf{A}$ and $\mathbf{B}$ can be met if $\min(n_1, n_2) = \sum_{i=1}^r L_i$. Those conditions may be met in unmixing applications however they are not ensured. We should note that the abovementioned uniqueness conditions were discussed for BTD without nonnegativity constraints. Nevertheless, the final solution always depends on the initialization due to the nonconvexity in blind unmixing. Matrix-Vector Nonnegative Tensor Factorization (MVNTF) was proposed in [105] which solves

$$ (\hat{\mathbf{K}}, \hat{\mathbf{B}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{K}, \mathbf{B}, \mathbf{E}} \frac{1}{2} \|\mathcal{Y} - \sum_{i=1}^{r} \mathbf{K}_i \mathbf{B}_i^T \odot \mathbf{e}_i\|_F^2 + \\ \frac{\lambda}{2} \| \sum_{i=1}^{r} \mathbf{K}_i \mathbf{B}_i^T - \mathbf{1}_{n_1} \mathbf{1}_{n_2}^T \|_F^2 \quad \text{s.t.} \quad \mathbf{K} \ge 0, \mathbf{B} \ge 0, \mathbf{E} \ge 0 \quad (45) $$

where the Frobenius norm for a third-order tensor is given by $\mathcal{X} = \sqrt{\sum_i \sum_j \sum_k x_{ijk}^2}$ and the second term ensures the ASC. Several variations of (45) were proposed in the literature. In [106], sparse and low-rank NTF was proposed which enforces both $\ell_1$ sparsity and nuclear norm on abundances. This method was further improved in [107] by exploiting weighted norms
for the abundances. In [108], the TV regularized MVNTF was proposed which utilizes the TV regularization on abundances to promote piecewise smoothness on the abundances and to capture spatial information. Tensor decomposition was also proposed for multi-feature hyperspectral unmixing to incorporate spatial features extracted by mathematical morphology [109].

4) Extended Linear Mixing Models and Spectral variability: An Extended Linear Mixing Model (ELMM) was proposed in [110] to address the endmember variability given by

$$ \mathbf{Y} = \mathbf{E}\mathbf{S}\mathbf{A} + \mathbf{N}, \mathbf{A} \geq 0, \mathbf{S} \geq 0, \mathbf{E} \geq 0 \quad (46) $$

where $\mathbf{S} = \mathbf{s}\mathbf{1}^T$ and $\mathbf{s}$ contains the scaling factor of the endmembers. Nonnegative constraint least square (NCLS) was used to estimate the scaled abundances. Assuming the scaling factor is the same for all the endmembers and knowing that the sum of the abundances is the scaling factor, therefore, the scaling factor was removed from abundances with a final normalization step. ELMM considers the spectral variation caused by illumination and topography. We should note that ELMM only considers uniform variations. In [111], a pixel-dependent ELMM was given by

$$ \mathbf{y}_i = \mathbf{E}_i \mathbf{a}_i + \mathbf{n}_i, \mathbf{A} \geq 0, \mathbf{E}_i = \mathbf{E} \odot \mathbf{c}_i \mathbf{1}^T, \mathbf{a}_i \geq 0, \mathbf{1}_r^T \mathbf{a}_i = 1, \quad (47) $$

where $\mathbf{C} \in \mathbb{R}^{r \times n}$ contains the scaling factors in its column. A minimization problem was proposed as

$$ \begin{aligned} (\hat{\mathbf{A}}, \hat{\mathcal{E}}, \hat{\mathcal{C}}) &= \arg \min_{\mathbf{A}, \mathcal{E}, \mathbf{C}} \frac{1}{2} \sum_{i=1}^{n} \|\mathbf{y}_i - \mathbf{E}_i \mathbf{a}_i\|_2^2 + \\ & \lambda_1 \sum_{i=1}^{n} \|\mathbf{E}_i - \mathbf{E} \odot \mathbf{c}_i \mathbf{1}^T\|_F^2 + \lambda_2 \psi(\mathbf{C}) + \lambda_3 \phi(\mathbf{A}) \quad (48) \\ \text{s.t.} \quad & \mathbf{A} \geq 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T \end{aligned} $$

where $\mathcal{E} \in \mathbb{R}^{p \times r \times n}$ is a $p \times r \times n$ tensor with $\mathbf{E}_i$ located in band $i$ and

$$ \psi(\mathbf{C}) = \frac{1}{2} (\|\mathbf{R}\mathbf{C}\|_F^2 + \|\mathbf{C}\mathbf{R}^T\|_F^2) \quad (49) $$

$$ \phi(\mathbf{A}) = \sum_{j=1}^{r} \|\mathbf{H}_h(\mathbf{a}_{(j)})\|_2 + \|\mathbf{H}_v(\mathbf{a}_{(j)})\|_2. \quad (50) $$

Generalized LMM (GLMM) was proposed in [112] to address the nonuniform variation. GLMM is given by

$$ \mathbf{y}_i = \mathbf{E}_i \mathbf{a}_i + \mathbf{n}_i, \text{ s.t. } \mathbf{E}_i = \mathbf{E} \odot \mathbf{C}_i, \mathbf{a}_i \geq 0, \mathbf{1}_r^T \mathbf{a}_i = 1, \quad (51) $$

where $\mathbf{C}_i \in \mathbb{R}^{p \times r}$, $i = 1, ..., n$ are scaling matrices dedicated to every pixel. The proposed minimization is given by

$$ \begin{aligned} (\hat{\mathbf{A}}, \hat{\mathcal{E}}, \hat{\mathcal{C}}) &= \arg \min_{\mathbf{A}, \mathcal{E}, \mathbf{C}} \frac{1}{2} \sum_{i=1}^{n} \|\mathbf{y}_i - \mathbf{E}_i \mathbf{a}_i\|_2^2 + \\ & \lambda_1 \sum_{i=1}^{n} \|\mathbf{E}_i - \mathbf{E} \odot \mathbf{C}_i\|_F^2 + \lambda_2 \sum_{i=1}^{n} \psi(\mathbf{C}_i) + \lambda_3 \phi(\mathbf{A}) \quad (52) \\ \text{s.t.} \quad & \mathbf{A} \geq 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \end{aligned} $$

where $\mathbf{C} \in \mathbb{R}^{p \times r \times n}$ is $p \times r \times n$ tensor with $\mathbf{C}_i$ located in band $i$. In [113], a framework was suggested to estimate the scaling matrices $\mathbf{C}_i$ using pure pixel endmember extraction

techniques. A perturbed linear mixing model (PLMM) was proposed in [114]. This simplifies the problem, however, the algorithm considerably relies on the presence of pure pixels and the extraction technique. The proposed model is given by

$$ \mathbf{y}_i = \mathbf{E}_i \mathbf{a}_i + \mathbf{n}_i, \text{ s.t. } \mathbf{E}_i = \mathbf{E} + \mathbf{C}_i \geq 0, \mathbf{E} \geq 0, \mathbf{a}_i \geq 0, \mathbf{1}_r^T \mathbf{a}_i = 1, \quad (53) $$

where $\mathbf{C}_i$ additive perturbation to model the spectral variability for every pixel. The proposed minimization is given by

$$ \begin{aligned} (\hat{\mathbf{A}}, \hat{\mathcal{E}}, \hat{\mathcal{C}}) &= \arg \min_{\mathbf{A}, \mathcal{E}, \mathbf{C}} \frac{1}{2} \sum_{i=1}^{n} \|\mathbf{y}_i - \mathbf{E}_i \mathbf{a}_i\|_2^2 + \frac{\lambda_1}{2} \|\mathbf{A}\mathbf{H}_d\|_F^2 \\ & \quad + \lambda_2 \psi(\mathbf{E}) + \lambda_3 \sum_{i=1}^{n} \|\mathbf{C}_i\|_F^2 \\ \text{s.t.} \quad & \mathbf{A} \geq 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \mathbf{E} \geq 0, \mathbf{E}_i \geq 0, i = 1, ..., n \quad (54) \end{aligned} $$

where $\psi$ can be $Vol^2$, $TV_e$, or the Euclidean distance between $\mathbf{E}$ and some reference spectral signatures and $\mathbf{H}_d$ is a matrix to calculate the deference of an abundance pixel with 4 nearest neighboring pixels (see [114] for more details). PLMM was used in [115], and GLMM was used in [116] for multitemporal unmixing. There exist other approaches that tackle spectral variability such as Sparsity-Enhanced Convolutional Decomposition (SeCoDe) [117] (a tensor-based model), Augmented Linear Mixing Models (ALMM) [118], and Subspace Unmixing with Low-Rank Attribute Embedding (SULoRA) [119].

## B. Shallow/Deep Neural Networks

Shallow or deep neural networks can be used for blind unmixing, as shown in Fig. 8. The decoder is often shallow, and a linear layer is often used for linear unmixing. In blind unmixing networks, the decoder weights are the network's parameters representing the endmembers to be learned. Therefore, the regularized loss function is given by

$$ \mathcal{L}_{reg}(\mathbf{y}, \hat{\mathbf{y}}) = \mathcal{L}(\mathbf{y}, \hat{\mathbf{y}}) + \lambda J(\mathbf{a}, \mathbf{W}_{en}, \mathbf{W}_{de}) \quad (55) $$

where $J$ can be a combination of several functions, such as a Frobenius norm on $\mathbf{W}_{en}$ with a sparsity norm such as $\ell_1$ norm on $\mathbf{a}$, and a Frobenius norm or a geometrical (MV induced) norm on $\mathbf{W}_{de}$. The loss function $\mathcal{L}$ might also be selected as a combination of several functions. For instance, EndNet [120] proposes a loss function with several terms, including a Kullback-Leibler divergence term, a SAD similarity, and a sparsity term, making the parameter selection very challenging.

The blind unmixing networks can be divided into two main groups: 1) autoencoder-based networks and 2) algorithm unrolling-based networks.

1) Autoencoder-based networks: The majority of the blind unmixing networks are variations of autoencoder-based networks. The main components of such networks were explained in Section III-D. However, the encoder can be much deeper and more complicated. Five main types of AEs were considered for blind unmixing [10]; 1) Sparse Nonnegative Autoencoders (SNAEs): SNAEs enforce sparsity on the bottleneck layer often via a regularization term in the loss, 2) Variational Autoencoders (VAEs): For VAEs, the decoder samples from a
for the abundances. In [108], the TV regularized MVNTF was proposed which utilizes the TV regularization on abundances to promote piecewise smoothness on the abundances and to capture spatial information. Tensor decomposition was also proposed for multi-feature hyperspectral unmixing to incorporate spatial features extracted by mathematical morphology [109].

4) Extended Linear Mixing Models and Spectral variability: An Extended Linear Mixing Model (ELMM) was proposed in [110] to address the endmember variability given by

$$ \mathbf{Y} = \mathbf{E}\mathbf{S}\mathbf{A} + \mathbf{N}, \mathbf{A} \geq 0, \mathbf{S} \geq 0, \mathbf{E} \geq 0 \quad (46) $$

where $\mathbf{S} = \mathbf{s}\mathbf{1}^T$ and $\mathbf{s}$ contains the scaling factor of the endmembers. Nonnegative constraint least square (NCLS) was used to estimate the scaled abundances. Assuming the scaling factor is the same for all the endmembers and knowing that the sum of the abundances is the scaling factor, therefore, the scaling factor was removed from abundances with a final normalization step. ELMM considers the spectral variation caused by illumination and topography. We should note that ELMM only considers uniform variations. In [111], a pixel-dependent ELMM was given by

$$ \mathbf{y}_i = \mathbf{E}_i\mathbf{a}_i + \mathbf{n}_i, \mathbf{A} \geq 0, \mathbf{E}_i = \mathbf{E} \odot \mathbf{c}_i\mathbf{1}^T, \mathbf{a}_i \geq 0, \mathbf{1}_r^T\mathbf{a}_i = 1, \quad (47) $$

where $\mathbf{C} \in \mathbb{R}^{r \times n}$ contains the scaling factors in its column. A minimization problem was proposed as

$$ \begin{aligned} (\hat{\mathbf{A}}, \hat{\mathcal{E}}, \hat{\mathcal{C}}) &= \arg \min_{\mathbf{A}, \mathcal{E}, \mathbf{C}} \frac{1}{2} \sum_{i=1}^{n} \|\mathbf{y}_i - \mathbf{E}_i \mathbf{a}_i\|_2^2 + \\ & \lambda_1 \sum_{i=1}^{n} \|\mathbf{E}_i - \mathbf{E} \odot \mathbf{c}_i \mathbf{1}^T\|_F^2 + \lambda_2 \psi(\mathbf{C}) + \lambda_3 \phi(\mathbf{A}) \quad (48) \\ \text{s.t.} \quad & \mathbf{A} \geq 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T \end{aligned} $$

where $\mathcal{E} \in \mathbb{R}^{p \times r \times n}$ is a $p \times r \times n$ tensor with $\mathbf{E}_i$ located in band $i$ and

$$ \psi(\mathbf{C}) = \frac{1}{2} (\|\mathbf{R}\mathbf{C}\|_F^2 + \|\mathbf{C}\mathbf{R}^T\|_F^2) \quad (49) $$

$$ \phi(\mathbf{A}) = \sum_{j=1}^{r} \|\mathbf{H}_h(\mathbf{a}_{(j)})\|_2 + \|\mathbf{H}_v(\mathbf{a}_{(j)})\|_2. \quad (50) $$

Generalized LMM (GLMM) was proposed in [112] to address the nonuniform variation. GLMM is given by

$$ \mathbf{y}_i = \mathbf{E}_i\mathbf{a}_i + \mathbf{n}_i, \text{ s.t. } \mathbf{E}_i = \mathbf{E} \odot \mathbf{C}_i, \mathbf{a}_i \geq 0, \mathbf{1}_r^T\mathbf{a}_i = 1, \quad (51) $$

where $\mathbf{C}_i \in \mathbb{R}^{p \times r}$, $i = 1, ..., n$ are scaling matrices dedicated to every pixel. The proposed minimization is given by

$$ \begin{aligned} (\hat{\mathbf{A}}, \hat{\mathcal{E}}, \hat{\mathcal{C}}) &= \arg \min_{\mathbf{A}, \mathcal{E}, \mathbf{C}} \frac{1}{2} \sum_{i=1}^{n} \|\mathbf{y}_i - \mathbf{E}_i \mathbf{a}_i\|_2^2 + \\ & \lambda_1 \sum_{i=1}^{n} \|\mathbf{E}_i - \mathbf{E} \odot \mathbf{C}_i\|_F^2 + \lambda_2 \sum_{i=1}^{n} \psi(\mathbf{C}_i) + \lambda_3 \phi(\mathbf{A}) \quad (52) \\ \text{s.t.} \quad & \mathbf{A} \geq 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \end{aligned} $$

where $\mathbf{C} \in \mathbb{R}^{p \times r \times n}$ is $p \times r \times n$ tensor with $\mathbf{C}_i$ located in band $i$. In [113], a framework was suggested to estimate the scaling matrices $\mathbf{C}_i$ using pure pixel endmember extraction

techniques. A perturbed linear mixing model (PLMM) was proposed in [114]. This simplifies the problem, however, the algorithm considerably relies on the presence of pure pixels and the extraction technique. The proposed model is given by

$$ \mathbf{y}_i = \mathbf{E}_i\mathbf{a}_i + \mathbf{n}_i, \text{ s.t. } \mathbf{E}_i = \mathbf{E} + \mathbf{C}_i \geq 0, \mathbf{E} \geq 0, \mathbf{a}_i \geq 0, \mathbf{1}_r^T\mathbf{a}_i = 1, \quad (53) $$

where $\mathbf{C}_i$ additive perturbation to model the spectral variability for every pixel. The proposed minimization is given by

$$ \begin{aligned} (\hat{\mathbf{A}}, \hat{\mathcal{E}}, \hat{\mathcal{C}}) &= \arg \min_{\mathbf{A}, \mathcal{E}, \mathbf{C}} \frac{1}{2} \sum_{i=1}^{n} \|\mathbf{y}_i - \mathbf{E}_i \mathbf{a}_i\|_2^2 + \frac{\lambda_1}{2} \|\mathbf{A}\mathbf{H}_d\|_F^2 \\ & \quad + \lambda_2 \psi(\mathbf{E}) + \lambda_3 \sum_{i=1}^{n} \|\mathbf{C}_i\|_F^2 \\ \text{s.t.} \quad & \mathbf{A} \geq 0, \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \mathbf{E} \geq 0, \mathbf{E}_i \geq 0, i = 1, ..., n \quad (54) \end{aligned} $$

where $\psi$ can be $Vol^2$, $TV_e$, or the Euclidean distance between $\mathbf{E}$ and some reference spectral signatures and $\mathbf{H}_d$ is a matrix to calculate the deference of an abundance pixel with 4 nearest neighboring pixels (see [114] for more details). PLMM was used in [115], and GLMM was used in [116] for multitemporal unmixing. There exist other approaches that tackle spectral variability such as Sparsity-Enhanced Convolutional Decomposition (SeCoDe) [117] (a tensor-based model), Augmented Linear Mixing Models (ALMM) [118], and Subspace Unmixing with Low-Rank Attribute Embedding (SULoRA) [119].

## B. Shallow/Deep Neural Networks

Shallow or deep neural networks can be used for blind unmixing, as shown in Fig. 8. The decoder is often shallow, and a linear layer is often used for linear unmixing. In blind unmixing networks, the decoder weights are the network's parameters representing the endmembers to be learned. Therefore, the regularized loss function is given by

$$ \mathcal{L}_{reg}(\mathbf{y}, \hat{\mathbf{y}}) = \mathcal{L}(\mathbf{y}, \hat{\mathbf{y}}) + \lambda J(\mathbf{a}, \mathbf{W}_{en}, \mathbf{W}_{de}) \quad (55) $$

where $J$ can be a combination of several functions, such as a Frobenius norm on $\mathbf{W}_{en}$ with a sparsity norm such as $\ell_1$ norm on $\mathbf{a}$, and a Frobenius norm or a geometrical (MV induced) norm on $\mathbf{W}_{de}$. The loss function $\mathcal{L}$ might also be selected as a combination of several functions. For instance, EndNet [120] proposes a loss function with several terms, including a Kullback-Leibler divergence term, a SAD similarity, and a sparsity term, making the parameter selection very challenging.

The blind unmixing networks can be divided into two main groups: 1) autoencoder-based networks and 2) algorithm unrolling-based networks.

1) Autoencoder-based networks: The majority of the blind unmixing networks are variations of autoencoder-based networks. The main components of such networks were explained in Section III-D. However, the encoder can be much deeper and more complicated. Five main types of AEs were considered for blind unmixing [10]; 1) Sparse Nonnegative Autoencoders (SNAEs): SNAEs enforce sparsity on the bottleneck layer often via a regularization term in the loss, 2) Variational Autoencoders (VAEs): For VAEs, the decoder samples from a
probability distribution model from the data and the encoded (hidden) features. The loss function uses a reconstructed loss plus a (KL-divergence) regularization term, ensuring that the encoded data stay close to a particular distribution (often a unit normal distribution), 3) Adversarial Autoencoders (AAEs): AAEs are based on generative adversarial networks (GANs). The network also uses a discriminator, and the AE's decoder often acts as a generator for the adversarial network. AAEs use a reconstructed plus an adversarial loss. 4) Denoising Autoencoders (DAEs): DAEs are trained based on partially corrupted input while minimizing the loss between the reconstructed and original (uncorrupted) inputs. Therefore, the network learns to remove the noise. 5) Convolutional Autoencoders (CAEs): CAEs use convolutional layers instead of fully connected layers. The weights are the convolutional filters that be learned throughout the training process. The convolutional filters capture the spatial correlation, which benefits the unmixing process. A comparative study was provided on the performance of AEs for hyperspectral unmixing in [10]. The probabilistic generative model for spectral unmixing (PGMSU) was proposed in [39]. The network uses a VAE to generate the endmembers and a parallel encoder to generate the abundances.

SNSA [121] utilizes a stack of nonnegative sparse autoencoders. The last one performs the task of unmixing, and the others are exploited to improve the robustness concerning the outliers. Deep AutoEncoder Network (DAEN) [122] exploits a stacked autoencoder to initialize a variational autoencoder which performs the unmixing task. An untied Denoising Autoencoder with Sparsity (uDAS) [123] exploits an additional denoising constraint on the decoder and a $\ell_{2,1}$ sparsity promoting constraint on the decoder. Note that the autoencoder receives a spectral pixel at a time for training the network. Training a network based on a single spectrum at a time ignores the spatial information. The advantage of incorporating spatial information for spectral unmixing has been confirmed in the literature. Therefore, patch-wise or cube-wise CNN was proposed to utilize the spatial information. In [124], an unmixing technique based on an autoencoder network is proposed that incorporates the spatial correlation between pixels by utilizing an adaptive abundance smoothing method. In [125], patch-wise or cube-wise convolutional autoencoders were proposed to incorporate the spatial information. In [126], parallel autoencoders were applied to spectral patches to exploit the spatial information. A CAE was proposed in [127]. The proposed network uses two convolutional layers as the decoder and a convolutional layer with linear activation as the decoder. ASC and ANC were enforced using a softmax, and SAD was used as the loss function. The cycle-consistency unmixing network [128] utilizes two convolutional autoencoders, which are cascaded and performed cyclically. The proposed loss function contains two terms for spectral reconstruction and one for abundance reconstruction to incorporate high-level semantic information. In [129], a combination of a CAE and a transformer was proposed for unmixing. The attention module of the transformer captures the nonlocal contextual feature dependencies between the image patches. In [130], VCA and FCLSU were first used to drive prior abundances

from superpixels. Then, the adversarial training guides the encoder of an AAE to transfer the spatial information into the network by matching the abundances. A VAE was used in [131] to address the spectral variability. First, a deep generative endmember model was learned using the pure pixel extracted from the data. Then, an optimization was suggested and solved using the endmember model as a prior. A deep generative model was proposed in [132], which uses three types of encoders, i.e., CNN, convolutional graph networks, and self-attention-based networks. A discriminator was used to enforce the generated endmembers to have a distribution similar to the endmember models either extracted from the data or selected from a library.

In AEs, a deeper encoder does not necessarily improve the performance of AE-based unmixing, however, it helps to enlarge the receptive field [10]. Therefore, in [37], a Multi-Stage CAE Network (MSNet) was suggested, which uses a multi-stage framework. Every stage uses a shallow encoder-decoder architecture. The input data are spatially downsampled at every stage. The features at the next stage were concatenated with the previous stage's input to preserve the abundance consistency of different stages.

A drawback of the many DL-based approaches is the absence of geometrical information. Most of them rely on SAD and/or MSE for training. As we already discussed and it has been proven [133], geometrical information could be very beneficial without pure pixels. To incorporate geometrical and spatial information, in [36], a minimum simplex convolutional network (MiSiCNet) was proposed. In [36], it was shown that the minimization problem

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \beta \|\mathbf{E} - \mathbf{m}\mathbf{1}_{r}^{\mathrm{T}}\|_{F}^{2} \\ + \lambda \phi(\mathbf{A}) \quad \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_{r}^{\mathrm{T}}\mathbf{A} = \mathbf{1}_{n}^{\mathrm{T}}, 0 \le \mathbf{E} \le 1 \tag{56} $$

can be turned into

$$ (\hat{\theta}_{1}, \hat{\theta}_{2}) = \arg \min_{\theta_{1}, \theta_{2}} \frac{1}{2} \|\mathbf{Y} - \theta_{2}f_{\theta_{1}}(\mathbf{Z})\|_{F}^{2} + \lambda \|\theta_{2} - \mathbf{m}\mathbf{1}_{r}^{\mathrm{T}}\|_{F}^{2} \\ \text{s.t.} \quad \hat{\mathbf{E}}\hat{\mathbf{A}} = \hat{\theta}_{2}f_{\hat{\theta}_{1}}(\mathbf{Z}) \tag{57} $$

where $\hat{\mathbf{Y}} = \hat{\mathbf{E}}\hat{\mathbf{A}}$, $\theta_1$ and $\theta_2$ are the parameters of the encoder and the decoder of the network, respectively. Therefore, the minimization problem (56) can be shifted to the minimization of the parameters of a deep network $f$ with a loss function given by:

$$ \mathcal{L}(\mathbf{Y}, \hat{\mathbf{Y}}, \hat{\mathbf{E}}, \mathbf{m}) = \frac{1}{2} \|\mathbf{Y} - \hat{\mathbf{Y}}\|_{F}^{2} + \lambda \|\hat{\mathbf{E}} - \mathbf{m}\mathbf{1}_{r}^{\mathrm{T}}\|_{F}^{2}, \tag{58} $$

where $f$ is an encoder-decoder network and $\hat{\mathbf{E}}$ is the weight of the decoder. MiSiCNet uses a softmax layer to enforce both the ASC and the ANC.

2) Algorithm Unrolling-based Networks: The main idea behind algorithm unrolling-based networks is to incorporate model-based with learning-based methods. In this group, an iterative algorithm (such as ADMM) is considered to solve a mixture model. Then, the iterations of the algorithm will be unrolled and implemented using neural networks. The updating equations will be substituted using NNs. Sequential
probability distribution model from the data and the encoded (hidden) features. The loss function uses a reconstructed loss plus a (KL-divergence) regularization term, ensuring that the encoded data stay close to a particular distribution (often a unit normal distribution), 3) Adversarial Autoencoders (AAEs): AAEs are based on generative adversarial networks (GANs). The network also uses a discriminator, and the AE's decoder often acts as a generator for the adversarial network. AAEs use a reconstructed plus an adversarial loss. 4) Denoising Autoencoders (DAEs): DAEs are trained based on partially corrupted input while minimizing the loss between the reconstructed and original (uncorrupted) inputs. Therefore, the network learns to remove the noise. 5) Convolutional Autoencoders (CAEs): CAEs use convolutional layers instead of fully connected layers. The weights are the convolutional filters that be learned throughout the training process. The convolutional filters capture the spatial correlation, which benefits the unmixing process. A comparative study was provided on the performance of AEs for hyperspectral unmixing in [10]. The probabilistic generative model for spectral unmixing (PGMSU) was proposed in [39]. The network uses a VAE to generate the endmembers and a parallel encoder to generate the abundances.

SNSA [121] utilizes a stack of nonnegative sparse autoencoders. The last one performs the task of unmixing, and the others are exploited to improve the robustness concerning the outliers. Deep AutoEncoder Network (DAEN) [122] exploits a stacked autoencoder to initialize a variational autoencoder which performs the unmixing task. An untied Denoising Autoencoder with Sparsity (uDAS) [123] exploits an additional denoising constraint on the decoder and a $\ell_{2,1}$ sparsity promoting constraint on the decoder. Note that the autoencoder receives a spectral pixel at a time for training the network. Training a network based on a single spectrum at a time ignores the spatial information. The advantage of incorporating spatial information for spectral unmixing has been confirmed in the literature. Therefore, patch-wise or cube-wise CNN was proposed to utilize the spatial information. In [124], an unmixing technique based on an autoencoder network is proposed that incorporates the spatial correlation between pixels by utilizing an adaptive abundance smoothing method. In [125], patch-wise or cube-wise convolutional autoencoders were proposed to incorporate the spatial information. In [126], parallel autoencoders were applied to spectral patches to exploit the spatial information. A CAE was proposed in [127]. The proposed network uses two convolutional layers as the decoder and a convolutional layer with linear activation as the decoder. ASC and ANC were enforced using a softmax, and SAD was used as the loss function. The cycle-consistency unmixing network [128] utilizes two convolutional autoencoders, which are cascaded and performed cyclically. The proposed loss function contains two terms for spectral reconstruction and one for abundance reconstruction to incorporate high-level semantic information. In [129], a combination of a CAE and a transformer was proposed for unmixing. The attention module of the transformer captures the nonlocal contextual feature dependencies between the image patches. In [130], VCA and FCLSU were first used to drive prior abundances

from superpixels. Then, the adversarial training guides the encoder of an AAE to transfer the spatial information into the network by matching the abundances. A VAE was used in [131] to address the spectral variability. First, a deep generative endmember model was learned using the pure pixel extracted from the data. Then, an optimization was suggested and solved using the endmember model as a prior. A deep generative model was proposed in [132], which uses three types of encoders, i.e., CNN, convolutional graph networks, and self-attention-based networks. A discriminator was used to enforce the generated endmembers to have a distribution similar to the endmember models either extracted from the data or selected from a library.

In AEs, a deeper encoder does not necessarily improve the performance of AE-based unmixing, however, it helps to enlarge the receptive field [10]. Therefore, in [37], a Multi-Stage CAE Network (MSNet) was suggested, which uses a multi-stage framework. Every stage uses a shallow encoder-decoder architecture. The input data are spatially downsampled at every stage. The features at the next stage were concatenated with the previous stage's input to preserve the abundance consistency of different stages.

A drawback of the many DL-based approaches is the absence of geometrical information. Most of them rely on SAD and/or MSE for training. As we already discussed and it has been proven [133], geometrical information could be very beneficial without pure pixels. To incorporate geometrical and spatial information, in [36], a minimum simplex convolutional network (MiSiCNet) was proposed. In [36], it was shown that the minimization problem

$$ (\hat{\mathbf{A}}, \hat{\mathbf{E}}) = \arg \min_{\mathbf{A}, \mathbf{E}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_{F}^{2} + \beta \|\mathbf{E} - \mathbf{m}\mathbf{1}_{r}^{\mathrm{T}}\|_{F}^{2} \\ + \lambda \phi(\mathbf{A}) \quad \text{s.t.} \quad \mathbf{A} \ge 0, \mathbf{1}_{r}^{\mathrm{T}}\mathbf{A} = \mathbf{1}_{n}^{\mathrm{T}}, 0 \le \mathbf{E} \le 1 \tag{56} $$

can be turned into

$$ (\hat{\theta}_{1}, \hat{\theta}_{2}) = \arg \min_{\theta_{1}, \theta_{2}} \frac{1}{2} \|\mathbf{Y} - \theta_{2}f_{\theta_{1}}(\mathbf{Z})\|_{F}^{2} + \lambda \|\theta_{2} - \mathbf{m}\mathbf{1}_{r}^{\mathrm{T}}\|_{F}^{2} \\ \text{s.t.} \quad \hat{\mathbf{E}}\hat{\mathbf{A}} = \hat{\theta}_{2}f_{\hat{\theta}_{1}}(\mathbf{Z}) \tag{57} $$

where $\hat{\mathbf{Y}} = \hat{\mathbf{E}}\hat{\mathbf{A}}$, $\theta_1$ and $\theta_2$ are the parameters of the encoder and the decoder of the network, respectively. Therefore, the minimization problem (56) can be shifted to the minimization of the parameters of a deep network $f$ with a loss function given by:

$$ \mathcal{L}(\mathbf{Y}, \hat{\mathbf{Y}}, \hat{\mathbf{E}}, \mathbf{m}) = \frac{1}{2} \|\mathbf{Y} - \hat{\mathbf{Y}}\|_{F}^{2} + \lambda \|\hat{\mathbf{E}} - \mathbf{m}\mathbf{1}_{r}^{\mathrm{T}}\|_{F}^{2}, \tag{58} $$

where $f$ is an encoder-decoder network and $\hat{\mathbf{E}}$ is the weight of the decoder. MiSiCNet uses a softmax layer to enforce both the ASC and the ANC.

2) Algorithm Unrolling-based Networks: The main idea behind algorithm unrolling-based networks is to incorporate model-based with learning-based methods. In this group, an iterative algorithm (such as ADMM) is considered to solve a mixture model. Then, the iterations of the algorithm will be unrolled and implemented using neural networks. The updating equations will be substituted using NNs. Sequential
repetitions of such networks will be trained using a loss function such as MSE. An advantage of such networks is that the tuning parameters can be trainable. SUnSAL with ADMM implementation was unrolled in [134] with an extra linear layer augmented to the network to implement blind unmixing (ADMMNet). The sequential networks in [134] were tested with and without parameter sharing. In [135], $ℓ_1$-NMF was considered for blind unmixing. The iterative soft-thresholding algorithm [136], was unrolled to estimate the abundances. An extra linear layer was used for the endmember estimation and added to the sequential networks. The scenarios of the networks with and without sharing parameters were also considered in [135]. A solution to the $ℓ_p$-NMF was unrolled for blind unmixing in [137]. Instead of using a linear layer, the endmember estimation step was also unrolled. Therefore, the sequential network was constructed of endmember networks followed by abundance networks. For the abundance and endmember networks, the generalized shrinkage thresholding solution given in [138] and the projected gradient descent were unrolled, respectively.

Finally, we should note that often blind NN-based unmixing can be turned into supervised unmixing by fixing the weights of the decoders using extracted endmembers.

### C. Statistical-based Approaches

Most recent advanced techniques developed for blind un-mixing fall into the first two categories. There are, however, other statistical-based approaches that cannot be reformulated as the constraint penalized least squares. The Bayesian ap-proach is widely used for blind unmixing. The maximum a posteriori (MAP) estimator and the minimum mean squared error (MMSE) estimator were used in [139] and [140], respec-tively. The Bayesian approaches using the MAP estimator can be reformulated as the penalized least square in the form of (30) when the random noise is assumed to be Gaussian and uncorrelated, i.e., the noise covariance matrix is diagonal [1], [141]. This assumption simplifies the inverse problem and its solution. Independent Component Analysis (ICA) [142], as a solution to blind source separation, was proposed for blind unmixing [143], [144]. The abundances are dependent (ASC), therefore, ICA is not a suitable solution to blind unmixing [145]. Dependent component analysis (DECA) was proposed for blind unmixing of highly mixed datasets [146], [147]. DECA assumes a mixture of Dirichlet densities as the prior for abundances. Dirichlet density satisfies ANC and ASC. The Bayesian estimation problems are often hard to solve, and the estimation of unknown parameters is computationally expensive. In [148], a Gibbs sampling strategy was proposed to decrease the complexity of the proposed hierarchical Bayesian model. A similar Bayesian approach was proposed in [33] for spectral mixture analysis (BayesianSMA) from Raman spectroscopy.

## VI. SEMISUPERVISED (LIBRARY-BASED) UNMIXING

When neither the pure pixels nor the sufficient spectra on the facets of the data simplex are available, the endmembers can not be successfully extracted/estimated, which results in

poor abundance estimations. Alternatively, blind unmixing techniques can be used. However, blind unmixing is a non-convex problem; therefore, a good initialization of endmem-bers using endmember extraction techniques will often benefit in finding a better solution. Additionally, in highly mixed scenarios, blind unmixing methods often fail due to the large sets of solutions that fit the data. Fig. 11 (d) demonstrates that in highly mixed scenarios, finding the real simplex is not trivial. Therefore, semi-supervised unmixing was suggested, relying on an endmember library. Here, we should note that we categorize all techniques that entirely or partially rely on a spectral library as semi-supervised unmixing techniques. Therefore, before discussing the semisupervised techniques, we discuss the library selection and construction.

Library selection or construction is a crucial step to the success of a semisupervised unmixing method. Blind selection of a library without extra attention and some processing steps will lead to poor results for semisupervised unmixing. There are two major paradigms in semisupervised unmixing i) Multiple Endmember Spectral and Mixture Analysis and ii) Sparse unmixing. The former was designed to address the spectral variability using the endmember variability. Therefore, the spectral library is designed to represent the variability of the endmembers. The latter seeks a sparse solution relying on the library, and therefore, the high correlation of the library endmembers avoids the sparse solution. In both paradigms, the library must well represent the materials in a sense, i.e., it must contain all the endmembers of the materials in the scene. Generally, a library can be obtained using either of the following ways [3], [149].

1) In situ field or/and laboratory measurements: endmem-ber libraries can be built using field or/and laboratory measurements [150]. Creating such libraries has several major drawbacks. It is hard, time-consuming, expensive, and might be sensor-dependent. Due to the different measuring conditions, systems, or instruments, there are often mismatches and scaling discrepancies between the endmembers from the library and the observed data [151], [152]. A list of available spectral libraries was given on https://specchio.ch/.

2) Construction using observed data: Extracting multiple pure pixels from the observed datasets representing every material in the observed dataset. These spectra can be further clustered to represent the endmember bundles. These libraries are data-dependent and often fail in highly mixed scenarios since there are not enough pure pixels in the dataset [61], [153], [154].

3) Construction using physical models: Radiative transfer models can be used to generate endmember libraries [155]. A radiative transfer model is defined as a function of physicochemical parameters. Therefore, different instances of endmembers can be obtained by varying physicochemical parameters. For instance, PROSPECT model [156] for vegetation and Hapke model [8] for reflectance modeling of intimate mixtures and densely packed grains or particles can be used. These libraries are sensor- and data-independent. However, they depend
repetitions of such networks will be trained using a loss function such as MSE. An advantage of such networks is that the tuning parameters can be trainable. SUnSAL with ADMM implementation was unrolled in [134] with an extra linear layer augmented to the network to implement blind unmixing (ADMMNet). The sequential networks in [134] were tested with and without parameter sharing. In [135], $\ell_1$-NMF was considered for blind unmixing. The iterative soft-thresholding algorithm [136], was unrolled to estimate the abundances. An extra linear layer was used for the endmember estimation and added to the sequential networks. The scenarios of the networks with and without sharing parameters were also considered in [135]. A solution to the $\ell_p$-NMF was unrolled for blind unmixing in [137]. Instead of using a linear layer, the endmember estimation step was also unrolled. Therefore, the sequential network was constructed of endmember networks followed by abundance networks. For the abundance and endmember networks, the generalized shrinkage thresholding solution given in [138] and the projected gradient descent were unrolled, respectively.

Finally, we should note that often blind NN-based unmixing can be turned into supervised unmixing by fixing the weights of the decoders using extracted endmembers.

### C. Statistical-based Approaches

Most recent advanced techniques developed for blind un-mixing fall into the first two categories. There are, however, other statistical-based approaches that cannot be reformulated as the constraint penalized least squares. The Bayesian ap-proach is widely used for blind unmixing. The maximum a posteriori (MAP) estimator and the minimum mean squared error (MMSE) estimator were used in [139] and [140], respec-tively. The Bayesian approaches using the MAP estimator can be reformulated as the penalized least square in the form of (30) when the random noise is assumed to be Gaussian and uncorrelated, i.e., the noise covariance matrix is diagonal [1], [141]. This assumption simplifies the inverse problem and its solution. Independent Component Analysis (ICA) [142], as a solution to blind source separation, was proposed for blind unmixing [143], [144]. The abundances are dependent (ASC), therefore, ICA is not a suitable solution to blind unmixing [145]. Dependent component analysis (DECA) was proposed for blind unmixing of highly mixed datasets [146], [147]. DECA assumes a mixture of Dirichlet densities as the prior for abundances. Dirichlet density satisfies ANC and ASC. The Bayesian estimation problems are often hard to solve, and the estimation of unknown parameters is computationally expensive. In [148], a Gibbs sampling strategy was proposed to decrease the complexity of the proposed hierarchical Bayesian model. A similar Bayesian approach was proposed in [33] for spectral mixture analysis (BayesianSMA) from Raman spectroscopy.

## VI. SEMISUPERVISED (LIBRARY-BASED) UNMIXING

When neither the pure pixels nor the sufficient spectra on the facets of the data simplex are available, the endmembers can not be successfully extracted/estimated, which results in

poor abundance estimations. Alternatively, blind unmixing techniques can be used. However, blind unmixing is a non-convex problem; therefore, a good initialization of endmem-bers using endmember extraction techniques will often benefit in finding a better solution. Additionally, in highly mixed scenarios, blind unmixing methods often fail due to the large sets of solutions that fit the data. Fig. 11 (d) demonstrates that in highly mixed scenarios, finding the real simplex is not trivial. Therefore, semi-supervised unmixing was suggested, relying on an endmember library. Here, we should note that we categorize all techniques that entirely or partially rely on a spectral library as semi-supervised unmixing techniques. Therefore, before discussing the semisupervised techniques, we discuss the library selection and construction.

Library selection or construction is a crucial step to the success of a semisupervised unmixing method. Blind selection of a library without extra attention and some processing steps will lead to poor results for semisupervised unmixing. There are two major paradigms in semisupervised unmixing i) Multiple Endmember Spectral and Mixture Analysis and ii) Sparse unmixing. The former was designed to address the spectral variability using the endmember variability. Therefore, the spectral library is designed to represent the variability of the endmembers. The latter seeks a sparse solution relying on the library, and therefore, the high correlation of the library endmembers avoids the sparse solution. In both paradigms, the library must well represent the materials in a sense, i.e., it must contain all the endmembers of the materials in the scene. Generally, a library can be obtained using either of the following ways [3], [149].

1) In situ field or/and laboratory measurements: endmem-ber libraries can be built using field or/and laboratory measurements [150]. Creating such libraries has several major drawbacks. It is hard, time-consuming, expensive, and might be sensor-dependent. Due to the different measuring conditions, systems, or instruments, there are often mismatches and scaling discrepancies between the endmembers from the library and the observed data [151], [152]. A list of available spectral libraries was given on https://specchio.ch/.

2) Construction using observed data: Extracting multiple pure pixels from the observed datasets representing every material in the observed dataset. These spectra can be further clustered to represent the endmember bundles. These libraries are data-dependent and often fail in highly mixed scenarios since there are not enough pure pixels in the dataset [61], [153], [154].

3) Construction using physical models: Radiative transfer models can be used to generate endmember libraries [155]. A radiative transfer model is defined as a function of physicochemical parameters. Therefore, different instances of endmembers can be obtained by varying physicochemical parameters. For instance, PROSPECT model [156] for vegetation and Hapke model [8] for reflectance modeling of intimate mixtures and densely packed grains or particles can be used. These libraries are sensor- and data-independent. However, they depend
on the availability and accuracy of a model for the target materials. Additionally, selecting representative samples from such complex models is a big challenge.

A combination of the libraries mentioned above also may be used. Alternatively, generative machine learning models can be used to augment a library. In [157], a variational autoencoder was used for library augmentation.

## A. Multiple Endmember Spectral and Mixture Analysis

As we discussed, endmember variability is a challenge in hyperspectral unmixing. A solution is to use a library containing endmembers' variations for materials. Then, the unmixing problem can be formulated as a semisupervised problem relying on such a dictionary. The pioneer method was proposed in [150] called multiple endmember spectral and mixture analysis (MESMA). MESMA was proposed to address the endmember variability. Assuming a structured library, $D = [D_1, D_2, ..., D_r]$ contains the endmember bundles, $D_i$'s, for all the materials, MESMA allows different and scaled endmembers for every pixel solving the problem

$$ \hat{\mathbf{a}}_i, \hat{\mathbf{E}}_i = \arg \min_{\mathbf{a}_i, \mathbf{E}_i} \frac{1}{2} \|\mathbf{y}_i - \mathbf{E}_i \mathbf{a}_i\|_2^2 \quad \text{s.t.} \\ \mathbf{E}_i \in D, \mathbf{a}_i \ge 0, \mathbf{1}_r^T \mathbf{a}_i = 1. \qquad (59) $$

This is a combinatorial problem, and MESMA searches for all combinations of endmembers (from every endmember bundle) and selects the one with the lowest reconstruction error. Endmember candidates are sometimes called models. This is solved for all the pixels. MESMA is a highly computationally demanding algorithm without any stopping criteria $\Pi_i m_i$ ($m_i$'s are dimensions of endmember bundles) times FCLSU must be computed for every pixel. However, this can be performed in parallel for all the pixels, it is computationally costly, and therefore, several algorithms have been proposed to reduce the computational burden. As a result, instead of FCLSU, the sum to one constrained least square unmixing (SCLSU) was used with a threshold for the error (RE). Note that ignoring the nonnegativity constraint leads to SCLSU, which has a closed-form solution and can be implemented efficiently. Then, solutions that are positive and less than one were selected [150], [158]. Additionally, the number of iterations and the threshold value can be changed to decrease the computations, but it can highly affect the performance [159]. In [160], Unconstrained Least Squares (UCLS) were used to speed up the search for the solutions to the combinatorial problem.

Alternatively, the size of the library can be decreased [161]–[163]. In [164], the angle between a pixel and the projections into hyperplanes of the fixed point (selected endmembers) and the endmember candidates were used to find the optimal endmember. The minimum angle is equivalent to the minimum error, and therefore, the SCLSU problem (note that projection into an affine set or a hyperplane only satisfies ASC) was shifted to more efficient angle computations. FCLSU was used to estimate the abundances using the selected endmembers (the optimal model), and the result with minimum reconstruction error was selected. On the other hand, in [157], a deep

generative model using variational autoencoders was suggested for augmenting the library endmembers to address the library mismatch. First, VAE was used to augment the library, and then MESMA was used for unmixing. However, this makes the algorithm computationally very expensive. A similar but more efficient approach was proposed in [165], which only uses one unmixing step.

## B. Sparse Unmixing

The first efficient idea was based on sparse and redundant modeling [31] and sparse regression, and therefore, it is often referred to as sparse unmixing. We should note that semisupervised unmixing refers to all the unmixing methods relying on a library and is not limited to sparse unmixing methods.

In sparse unmixing, the fractional abundances are estimated using sparse regression techniques. These methods describe each spectrum as a sparse linear combination of the elements of a rich library of pure spectra.

The main idea of sparse unmixing was first suggested in [31]. If a well-designed dictionary is available, the sparse regression can be used to estimate the abundances given by

$$ \hat{\mathbf{X}} = \arg \min_{\mathbf{X}} \frac{1}{2} \|\mathbf{Y} - \mathbf{D}\mathbf{X}\|_F^2 + \lambda \|\mathbf{X}\|_1 \\ \text{s.t.} \quad \mathbf{X} \ge 0, \mathbf{1}_m^T \mathbf{X} = \mathbf{1}_n^T. \qquad (60) $$

An ADMM-based algorithm was proposed to solve the problem (60), and therefore it was called sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL). In [31], also the constraint SUnSAL (CSUnSAL) proposed given by

$$ \hat{\mathbf{X}} = \arg \min_{\mathbf{X}} \|\mathbf{X}\|_1 \quad \text{s.t.} \quad \|\mathbf{Y} - \mathbf{D}\mathbf{X}\|_F < \delta, \\ \mathbf{X} \ge 0, \mathbf{1}_m^T \mathbf{X} = \mathbf{1}_n^T, \qquad (61) $$

and was solved using an ADMM-based algorithm. However, it is suggested to use SUnSAL without ASC due to the conflict with $\ell_1$ [90]. Additionally, ASC was found to be a rigorous constraint that often does not occur in the real world due to noise and signature variability [152]. Therefore, SUnSAL often refers to the problem (60) without ASC, and problem (60) refers to FCLS version of SUnSAL which can be seen in the arXiv version of SUnSAL [166]. In [167], problem (60) is suggested for implementation on a GPU when $\lambda = 0$. In [152], the nonconvex $\ell_0$ sparse regression

$$ \hat{\mathbf{X}} = \arg \min_{\mathbf{X}} \|\mathbf{X}\|_0 \quad \text{s.t.} \quad \|\mathbf{Y} - \mathbf{D}\mathbf{X}\|_F < \delta, \\ \mathbf{X} \ge 0, \qquad (62) $$

were compared with SUnSAL and CSUnSAL. The nonnegative $\ell_0$ sparse regression problem (62) was solved using a greedy algorithm called orthogonal matching pursuit [168]. It was shown that SUnSAL outperforms the other two sparse regressions. Additionally, the experiments in [152] confirmed the advantage of ANC.

The success of sparse unmixing depends on the following assumptions; i) the library must well represent the materials in
on the availability and accuracy of a model for the target materials. Additionally, selecting representative samples from such complex models is a big challenge.

A combination of the libraries mentioned above also may be used. Alternatively, generative machine learning models can be used to augment a library. In [157], a variational autoencoder was used for library augmentation.

A. Multiple Endmember Spectral and Mixture Analysis

As we discussed, endmember variability is a challenge in hyperspectral unmixing. A solution is to use a library containing endmembers' variations for materials. Then, the unmixing problem can be formulated as a semisupervised problem relying on such a dictionary. The pioneer method was proposed in [150] called multiple endmember spectral and mixture analysis (MESMA). MESMA was proposed to address the endmember variability. Assuming a structured library, $D = [D_1, D_2, ..., D_r]$ contains the endmember bundles, $D_i$'s, for all the materials, MESMA allows different and scaled endmembers for every pixel solving the problem

$$ 
\begin{gathered}
\hat{\mathbf{a}}_i, \hat{\mathbf{E}}_i = \arg \min_{\mathbf{a}_i, \mathbf{E}_i} \frac{1}{2} \|\mathbf{y}_i - \mathbf{E}_i \mathbf{a}_i\|_2^2 \quad \text{s.t.} \\
\mathbf{E}_i \in D, \mathbf{a}_i \ge 0, \mathbf{1}_r^T \mathbf{a}_i = 1.
\end{gathered}
\quad (59) $$

This is a combinatorial problem, and MESMA searches for all combinations of endmembers (from every endmember bundle) and selects the one with the lowest reconstruction error. Endmember candidates are sometimes called models. This is solved for all the pixels. MESMA is a highly computationally demanding algorithm without any stopping criteria $\Pi_i m_i$ ($m_i$'s are dimensions of endmember bundles) times FCLSU must be computed for every pixel. However, this can be performed in parallel for all the pixels, it is computationally costly, and therefore, several algorithms have been proposed to reduce the computational burden. As a result, instead of FCLSU, the sum to one constrained least square unmixing (SCLSU) was used with a threshold for the error (RE). Note that ignoring the nonnegativity constraint leads to SCLSU, which has a closed-form solution and can be implemented efficiently. Then, solutions that are positive and less than one were selected [150], [158]. Additionally, the number of iterations and the threshold value can be changed to decrease the computations, but it can highly affect the performance [159]. In [160], Unconstrained Least Squares (UCLS) were used to speed up the search for the solutions to the combinatorial problem.

Alternatively, the size of the library can be decreased [161]–[163]. In [164], the angle between a pixel and the projections into hyperplanes of the fixed point (selected endmembers) and the endmember candidates were used to find the optimal endmember. The minimum angle is equivalent to the minimum error, and therefore, the SCLSU problem (note that projection into an affine set or a hyperplane only satisfies ASC) was shifted to more efficient angle computations. FCLSU was used to estimate the abundances using the selected endmembers (the optimal model), and the result with minimum reconstruction error was selected. On the other hand, in [157], a deep

generative model using variational autoencoders was suggested for augmenting the library endmembers to address the library mismatch. First, VAE was used to augment the library, and then MESMA was used for unmixing. However, this makes the algorithm computationally very expensive. A similar but more efficient approach was proposed in [165], which only uses one unmixing step.

B. Sparse Unmixing

The first efficient idea was based on sparse and redundant modeling [31] and sparse regression, and therefore, it is often referred to as sparse unmixing. We should note that semisupervised unmixing refers to all the unmixing methods relying on a library and is not limited to sparse unmixing methods.

In sparse unmixing, the fractional abundances are estimated using sparse regression techniques. These methods describe each spectrum as a sparse linear combination of the elements of a rich library of pure spectra.

The main idea of sparse unmixing was first suggested in [31]. If a well-designed dictionary is available, the sparse regression can be used to estimate the abundances given by

$$ 
\begin{gathered}
\hat{\mathbf{X}} = \arg \min_{\mathbf{X}} \frac{1}{2} \|\mathbf{Y} - \mathbf{D}\mathbf{X}\|_F^2 + \lambda \|\mathbf{X}\|_1 \\
\text{s.t.} \quad \mathbf{X} \ge 0, \mathbf{1}_m^T \mathbf{X} = \mathbf{1}_n^T.
\end{gathered}
\quad (60) $$

An ADMM-based algorithm was proposed to solve the problem (60), and therefore it was called sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL). In [31], also the constraint SUnSAL (CSUnSAL) proposed given by

$$ 
\begin{gathered}
\hat{\mathbf{X}} = \arg \min_{\mathbf{X}} \|\mathbf{X}\|_1 \quad \text{s.t.} \quad \|\mathbf{Y} - \mathbf{D}\mathbf{X}\|_F < \delta, \\
\mathbf{X} \ge 0, \mathbf{1}_m^T \mathbf{X} = \mathbf{1}_n^T,
\end{gathered}
\quad (61) $$

and was solved using an ADMM-based algorithm. However, it is suggested to use SUnSAL without ASC due to the conflict with $\ell_1$ [90]. Additionally, ASC was found to be a rigorous constraint that often does not occur in the real world due to noise and signature variability [152]. Therefore, SUnSAL often refers to the problem (60) without ASC, and problem (60) refers to FCLS version of SUnSAL which can be seen in the arXiv version of SUnSAL [166]. In [167], problem (60) is suggested for implementation on a GPU when $\lambda = 0$. In [152], the nonconvex $\ell_0$ sparse regression

$$ 
\begin{gathered}
\hat{\mathbf{X}} = \arg \min_{\mathbf{X}} \|\mathbf{X}\|_0 \quad \text{s.t.} \quad \|\mathbf{Y} - \mathbf{D}\mathbf{X}\|_F < \delta, \\
\mathbf{X} \ge 0,
\end{gathered}
\quad (62) $$

were compared with SUnSAL and CSUnSAL. The nonnegative $\ell_0$ sparse regression problem (62) was solved using a greedy algorithm called orthogonal matching pursuit [168]. It was shown that SUnSAL outperforms the other two sparse regressions. Additionally, the experiments in [152] confirmed the advantage of ANC.

The success of sparse unmixing depends on the following assumptions; i) the library must well represent the materials in