# Fast Semisupervised Unmixing Using Nonconvex Optimization

Behnood Rasti, Alexandre Zouaoui, Julien Mairal, Jocelyn Chanussot

► To cite this version:

Behnood Rasti, Alexandre Zouaoui, Julien Mairal, Jocelyn Chanussot. Fast Semisupervised Unmixing Using Nonconvex Optimization. IEEE Transactions on Geoscience and Remote Sensing, 2024, 62, pp.5526713:1-13. 10.1109/TGRS.2024.3440663. hal-04409409v2

HAL Id: hal-04409409

https://hal.science/hal-04409409v2

Submitted on 27 Sep 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
# Fast Semi-supervised Unmixing using Non-convex Optimization

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

***Abstract***—In this paper, we introduce a novel linear model tailored for semisupervised/library-based unmixing. Our model incorporates considerations for library mismatch while enabling the enforcement of the abundance sum-to-one constraint (ASC). Unlike conventional sparse unmixing methods, this model involves nonconvex optimization, presenting significant computational challenges. We demonstrate the efficacy of Alternating Methods of Multipliers (ADMM) in cyclically solving these intricate problems. We propose two semisupervised unmixing approaches, each relying on distinct priors applied to the new model in addition to the ASC: sparsity prior and convexity constraint. Our experimental results validate that enforcing the convexity constraint outperforms the sparsity prior for the endmember library. These results are corroborated across three simulated datasets (accounting for spectral variability and varying pixel purity levels) and the Cuprite dataset. Additionally, our comparison with conventional sparse unmixing methods showcases considerable advantages of our proposed model, which entails nonconvex optimization. Notably, our implementations of the proposed algorithms—fast semisupervised unmixing (FaSUN) and sparse unmixing using soft-shrinkage (SUNs)—prove considerably more efficient than traditional sparse unmixing methods. SUNS and FaSUN were implemented using PyTorch and provided in a dedicated Python package called Fast Semisupervised Unmixing (FUNmix), which is open-source and available at https://github.com/BehnoodRasti/FUnmix.

***Index Terms***—sparse unmixing, hyperspectral, sparsity, semi-supervised, blind, unmixing, PyTorch, GPU, alternating direction method of multipliers, nonconvex, optimization

## I. INTRODUCTION

SPECTRAL unmixing estimates the abundances of pure spectra of materials called endmembers. Depending on the prior knowledge of endmembers, the unmixing problem can be categorized into three main types: supervised unmixing, blind unmixing, and semi-supervised unmixing. In supervised unmixing, abundances are estimated relying on known endmembers. Blind unmixing estimates both endmembers and abundances simultaneously. Semi-supervised unmixing relies on an endmember library to estimate the corresponding abundances [1].

An unmixing problem can be tackled using a sequential process of extracting/estimating endmembers (often using a

Behnood Rasti (corresponding author) is with the Faculty of Electrical Engineering and Computer Science, Technische Universität Berlin, and the Berlin Institute for the Foundations of Learning and Data, Berlin, Germany, behnood.rasti@tu-berlin.de, 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 2024; revised 2024.

geometrical approach) followed by an abundance estimation technique. We call this group of methods supervised since the endmembers are assumed to be known in the abundance estimating step [1]. Endmember extraction/estimation techniques often rely on the geometry of the data simplex, such as Vertex Component Analysis (VCA) [2], Simplex volume maximization (SiVM) [3], the minimum volume simplex analysis (MVSA) [4], and the simplex identification via variable splitting and augmented Lagrangian (SISAL) [5]. The discussion on the endmember extraction/estimation techniques is out of the scope of this paper, and therefore, we refer to [6], [1] for an overview of this topic.

When the endmembers are known, abundances can be estimated. Unconstrained least squares unmixing (UCLSU) via the orthogonal subspace projection was proposed for abundance estimation [7]. Non-negative constrained least squares unmixing (NCLSU) [8], [9] was proposed to estimate the abundances subjected to ANC. There are several attempts to solve the least squares problem subjected to both ANC and ASC [10], [11]. The first efficient algorithm was proposed in [12] and called fully constrained least squares unmixing (FCLSU). FCLSU can be efficiently solved using general-purpose convex optimization toolboxes. However, in this paper, we will show that with the advances in graphical processing units (GPU), FCLSU can be efficiently solved using convex optimization techniques i.e., Alternating Direction Method of Multipliers (ADMM) [13].

The pioneer semisupervised unmixing is Multiple Endmember Spectral and Mixture Analysis (MESMA) [14] proposed to address endmember variability. MESMA assumes a structured library containing endmember bundles for all materials, allowing different scaled endmembers for each pixel. However, this is a combinatorially complex task and computationally expensive. J. M. Bioucas Dias and M. A. T. Figueiredo [15] proposed a formulation of the semi-supervised unmixing problem as a sparse regression problem, giving rise to what is known as sparse unmixing. They proposed sparse unmixing by variable splitting and augmented lagrangian (SUNSAL), as well as its variant known as Constrained SUNSAL (C-SUNSAL) [15]. SUNSAL and C-SUNSAL employ the $l_1$ penalty to encourage sparsity in the abundance estimation. In the case of SUNSAL, the $l_2$ norm is combined with the $l_1$ penalty to enhance fidelity, whereas C-SUNSAL uses the $l_2$ norm as a constraint in order to minimize the $l_1$ term. The optimization problems associated with SUNSAL were tackled through the Alternating Direction Method of Multipliers (ADMM) [13].

Sparse unmixing offers computational efficiency; however,
the high correlation among library endmembers presents a significant challenge for sparse regression. This concern has been addressed through techniques like library pruning and the application of various sparsifying regularizers. To enhance SUnSAL, a total variation (TV) penalty was incorporated, resulting in SUnSAL-TV [16], which leverages spatial information. Nevertheless, SUnSAL-TV does not satisfy the constraint of abundance sum-to-one (ASC) due to its potential conflict with the $l_1$ penalty.

Collaborative sparse unmixing [17] enforces a constraint by applying the sum of $l_2$ norms to the abundances. Double Reweighted Sparse Unmixing (DRSU) [18] and Spectral-Spatial Weighted Sparse Unmixing (S2WSU) [19] adopt a weighted $l_1$ norm approach to induce sparsity. Additionally, DRSU employs total variation to capture spatial attributes. The Multiscale Sparse Unmixing Algorithm (MUA) [20] capitalizes on spatial correlation by performing sparse regression on segmented pixels, enabling the capture of both spectral variability and spatial correlation. In this context, segmentation techniques such as Binary Partition Tree (BPT), Simple Linear Iterative Clustering (SLIC), and the K-means algorithm were suggested in [20]. In [21], SLIC was adopted for segmentation, while sparse unmixing was executed with superpixel-based graph Laplacian regularization.

A common drawback of the aforementioned sparse unmixing techniques is that the estimated fractional abundances may not accurately represent the aerial fraction of each pure material on the ground due to the absence of the ASC constraint. As mentioned, applying $l_1$ penalties to the abundances cannot maintain the ASC. This issue was addressed in sparse unmixing using a convolutional neural network (SUnCNN) [22]. In [22], we demonstrated that selecting an appropriate prior for sparse regression could be transformed into an optimization task involving the parameters of a deep encoder-decoder network, while the ASC could be enforced using a softmax layer. However, it is worth noting that selecting suitable hyperparameters for such a deep network is often a challenging endeavor. In [23], an asymmetric encoder-decoder network is used with a sparse variation of softmax to avoid the full support of softmax while enforcing ASC.

Algorithm unrolling-based strategies have also been explored in the context of sparse unmixing. A recent instance of this approach involves the development of a shallow network for sparse unmixing, as outlined in [24] and [25]. In these works, an unrolling technique was employed to address the nonnegative $l_1$ sparse regression problem, i.e., SUnSAL. To enhance spatial information integration, an intermediate convolutional layer was applied to the abundance representation. The training of the shallow network involved a combination of loss functions, including SAD (Sum of Absolute Differences), MSE (Mean Squared Error), and SID (Spectral Information Divergence).

Furthermore, in [26], a similar unrolling approach was employed, where the Iterative Soft-Thresholding Algorithm (ISTA) [27] was unrolled to tackle the nonnegative $l_1$ sparse regression problem. Additionally, in the pursuit of sparse unmixing, the unrolling technique was applied to SUnSAL, as documented in [28].

A significant concern arises with the techniques mentioned above when the endmembers do not align with those in the library. In such cases, even a well-curated and pruned spectral library may fall short in representing all the unique endmembers present in real-world datasets. Factors such as noise, atmospheric effects, variations in illumination, and intrinsic material differences introduce shifts in the endmembers. This often leads to scaling discrepancies between the endmembers in the scene and those in the library.

We recently introduced Sparse Unmixing using Archetypal Analysis (SUnAA) to tackle this issue. SUnAA assumes that the endmembers can be expressed as convex combinations of the library endmembers. It addresses the problem as a non-convex optimization using a cyclic descent algorithm. SUnAA runs on the CPU and can be computationally demanding. The number of pixels and endmembers in the dataset influences its performance. In this paper, we propose efficient unmixing algorithms using ADMM to address those drawbacks. Overall, our main contributions can be summarized in three key aspects:

1) We introduce a new linear model for semisupervised/library-based unmixing which takes into account the endmember library mismatch and ASC. Our experiments show the advantages of this model compared with the sparse and redundant model used in conventional sparse unmixing.

2) We propose two ADMM-based methods named Fast Semisupervised Unmixing (FaSUn) and Fast Sparse Unmixing using Soft-Shrinkage (SUnS) aiming at comparing two different priors on the new model. FaSUn and SUnS enforce convexity and sparsity on the endmembers, respectively. Our experiments reveal that the convexity constraint outperforms the sparsity prior.

3) We provide GPU (PyTorch)-based implementations for the ADMM-based algorithms showcasing the efficiency of the proposed algorithms compared to the state-of-the-art semisupervised unmixing techniques.

## II. METHODOLOGY

### A. Low-rank Linear Mixture Model

Assuming matrix $E \in \mathbb{R}^{p \times r}$ contains $r$ endmembers within the observed hyperspectral pixel $y \in \mathbb{R}^p$ (i.e., the sensor has $p$ bands), then linear mixture model (LMM) is given by

$$ 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. To represent all the pixels we use the matrix notation $Y$. Then, we have [29]

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

where $Y \in \mathbb{R}^{p \times n}$ is the observed HSI, with $n$ pixels and $p$ bands, $N \in \mathbb{R}^{p \times n}$ is noise, and $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. LMM is often used for supervised modeling [1].
## B. Sparse and Redundant Linear Mixture Model

Sparse and Redundant Linear Mixture Model is given by

$$ 
\begin{gathered}
\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,
\end{gathered}
\quad (3) $$

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. Please note that $\mathbf{D}$ serves as an overcomplete dictionary, and as such, it should be meticulously crafted. A well-structured dictionary comprises endmembers representing the materials present in the scene and can efficiently reduce the redundancy in $\mathbf{X}$. Consequently, it becomes possible to prune a spectral library based on the spectral angles between spectra, meaning that spectra with small angular differences are removed. However, there is a caveat: this pruning strategy carries the risk of losing endmember materials if they happen to be scaled versions of each other.

In the context of a well-designed dictionary, the pixels in the scene are composed of a mixture of a few dictionary atoms. This characteristic results in $\mathbf{X}$ being a sparse matrix. It is worth noting that if a specific endmember material is absent from the observed spectra, the corresponding row in $\mathbf{X}$ will be entirely composed of zeros. This is a common occurrence since abundance values are typically sparse in this model. This framework is frequently employed in the context of sparse unmixing, where fractional abundances $\mathbf{X}$ are estimated by applying sparsity-enforcing penalties or constraints within a sparse regression formulation.

## C. A New Linear Mixture Model

To enforce both ASC and sparsity, recently, a new linear model inspired by archetypal analysis [30] was proposed in [29] for library-based unmixing (semisupervised). In [29], we proposed a mixing model for observed spectra in

$$ \mathbf{Y} = \mathbf{D}\mathbf{B}\mathbf{A} + \mathbf{N}, \quad (4) $$

$\mathbf{B} \in \mathbb{R}^{m \times r}$, determines the contributions of the endmembers from $\mathbf{D}$. Model (4) exploits both low-rank property and sparse contribution of the endmembers. The low-rank property of model (4) decreases the computational time. However, the problem will turn into a non-convex problem (both $\mathbf{B}$ and $\mathbf{A}$ are unknown), increasing the algorithms' complexity. There are two main advantages of the model (4) compared to the sparse and redundant model: 1- The ASC can be enforced, 2- It can compensate for the mismatch between the library endmembers and data endmembers [29]. Incorporating matrix $\mathbf{B}$ into the model formulation allows the model to account for the mismatch between data endmembers and library endmembers. The classical sparse and redundant LMM, i.e., $\mathbf{D}\mathbf{X}$, assumes that the library endmembers should match the data endmembers. However, in a real application, the actual endmembers that exist in the scene differ from the library ones (endmember variability). Applying different prior on matrix $\mathbf{B}$ allows different linear combinations of the library endmembers to represent the data endmembers. On the other hand, in model

(4) the number of endmembers is unknown and challenging to estimate. 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 [31], [32], [33]. In practice, this problem can be addressed using eigenvalue-based detection techniques [34], [33], [35] or estimating the mean square errors [31], [32]. Geometrical-based approaches were also proposed for endmember estimation [36], [37], [38]. Indeed, this is an advantage of the conventional sparse and redundant model (3).

## D. FaSUn: Fast Semisupervised Unmixing

We propose a nonconvex optimization to simultaneously estimate $\mathbf{B}$ and $\mathbf{A}$:

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

Note that, in (5), the unknown endmembers are a convex combination of the library's endmembers due to the non-negativity and sum to one constraint on $\mathbf{B}$. In [29], we proposed a parameter-free solution to (5) using active set methods. The major issue with SUnAA is that it runs on the CPU and is highly time-consuming. Here, we propose an ADMM-based solution for the proposed minimization problem (5), which benefits from a GPU-accelerated implementation.

First, we should note that the minimization problem (5) is non-convex; however, it can be solved in two steps using a cyclic descent scheme; A-step: when $\mathbf{B}$ is fixed and B-step: $\mathbf{A}$ is fixed. In every step, we are dealing with a convex optimization, and therefore, every solution of the steps successively decreases the loss function, which leads to a minimum. The convergence of the final solution is guaranteed upon the convergence of every step throughout the iterations.

A-step: when $\mathbf{B}$ is fixed then $\mathbf{E} = \mathbf{D}\mathbf{B}$ is fixed. Therefore, problem 5 turns to

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

Problem (6) can be solved using any convex optimization, least squares, or quadratic programming solver. However, for unmixing and particularly Earth observation applications, we are often dealing with big datasets, and therefore, these general-purpose convex optimization solvers are not efficient. Here, we propose an ADMM solution.

To solve problem (6), we start by splitting $\mathbf{A}$,

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

Using ADMM, the augmented Lagrangian (AL) can be written as

$$ 
\begin{gathered}
\hat{\mathbf{A}}, \hat{\mathbf{S}} = \arg \min_{\mathbf{A},\mathbf{S}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \frac{\mu}{2} \|\mathbf{S} - \mathbf{A} - \mathbf{L}\|_F^2 \\
\text{s.t.} \quad \mathbf{S} \ge 0, \mathbf{1}_r^T\mathbf{A} = \mathbf{1}_n^T,
\end{gathered}
\quad (8) $$
where L is the Lagrange multiplier. Note that, we did not use AL for ASC. As can be seen in Appendix A, using Lagrangian and Karush-Kuhn-Tucker (KKT) will lead to a close form solution and the augmented term in (8) turns the matrix needs to be inverted non-singular. The solution to this problem is given in three steps.

When S is fixed, the problem turns to

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

Problem (9) is a quadratic programming (also known as least squares) with equality constraint (QuEC). In appendix A we show that there is a closed-form solution for (9) which is given by

$$ \mathbf{A} = \mathrm{QuEC}(\mathbf{A}, \mathbf{S}, \mathbf{L}; \mathbf{Y}, \mathbf{E}, \mu), \quad (10) $$

where QuEC is the function given by

$$ \mathrm{QuEC}(\mathbf{A}, \mathbf{S}, \mathbf{L}; \mathbf{Y}, \mathbf{E}, \mu) = (\mathbf{Q} + \mathbf{Q}_1 \mathbf{r} \mathbf{c}_1 \mathbf{r}^T \mathbf{Q})(\mathbf{E}^T \mathbf{Y} + \mu(\mathbf{S} - \mathbf{L})) - \mathbf{Q}_1 \mathbf{r} \mathbf{c}_1 \mathbf{n}^T, \quad (11) $$

and

$$ 
\begin{aligned}
\mathbf{Q} &= (\mathbf{E}^T \mathbf{E} + \mu \mathbf{I}_r)^{-1} \\
c &= -1/(\mathbf{1}_r^T \mathbf{Q} \mathbf{1}_r),
\end{aligned}
\quad (12) $$

When A is fixed, the problem turns to

$$ \hat{\mathbf{S}} = \arg \min_{\mathbf{S}} \frac{\mu}{2} \|\mathbf{S} - \mathbf{A} - \mathbf{L}\|_F^2 \quad \text{s.t.} \quad \mathbf{S} \ge 0, $$

and the solution is given by

$$ \mathbf{S} = \max(0, \mathbf{A} + \mathbf{L}), \quad (13) $$

Finally, we update the multiplier

$$ \mathbf{L} = \mathbf{L} + \mathbf{A} - \mathbf{S}. \quad (14) $$

B-step: When A is fixed, problem (5) turns to

$$ 
\begin{aligned}
\hat{\mathbf{B}} &= \arg \min_{\mathbf{B}} \frac{1}{2} \|\mathbf{Y} - \mathbf{D}\mathbf{B}\mathbf{A}\|_F^2 \\
\text{s.t.} \quad \mathbf{B} \ge 0, \quad \mathbf{1}_m^T \mathbf{B} = \mathbf{1}_r^T.
\end{aligned}
\quad (15) $$

Splitting the variables as $\mathbf{B} = \mathbf{S}_1$ and $\mathbf{D}\mathbf{B} = \mathbf{S}_2$, the AL is given by

$$ 
\begin{aligned}
\hat{\mathbf{B}} = & \arg \min_{\mathbf{B}, \mathbf{S}_1, \mathbf{S}_2} \frac{1}{2} \|\mathbf{Y} - \mathbf{S}_2 \mathbf{A}\|_F^2 + \frac{\mu}{2} \|\mathbf{S}_1 - \mathbf{B} - \mathbf{L}_1\|_F^2 \\
& + \frac{\mu_1}{2} \|\mathbf{S}_2 - \mathbf{D}\mathbf{B} - \mathbf{L}_2\|_F^2 \quad \text{s.t.} \quad \mathbf{S}_2 \ge 0, \quad \mathbf{1}_m^T \mathbf{B} = \mathbf{1}_r^T.
\end{aligned}
\quad (16) $$

We solve (16) with respect to each unknown matrix separately. Therefore, we have

$$ \hat{\mathbf{B}} = \mathrm{QuEC}(\mathbf{S}_1, \mathbf{L}_1; (\mathbf{S}_2 - \mathbf{L}_2), \mathbf{D}, \mu_1/\mu_2), \quad (17) $$

$$ \hat{\mathbf{S}}_1 = \max(0, \mathbf{B} + \mathbf{L}_1), \quad (18) $$

$$ \hat{\mathbf{S}}_2 = (\mathbf{Y}\mathbf{A}^T + \mu_2(\mathbf{D}\mathbf{B} + \mathbf{L}_2))(\mathbf{A}\mathbf{A}^T + \mu_2\mathbf{I}_r)^{-1}. \quad (19) $$

Finally, we update the multipliers

$$ \mathbf{L}_1 = \mathbf{L}_1 + \mathbf{B} - \mathbf{S}_1, \quad (20) $$

$$ \mathbf{L}_2 = \mathbf{L}_2 + \mathbf{D}\mathbf{B} - \mathbf{S}_2. \quad (21) $$

### Algorithm 1: FASUn

**Input:** $\mathbf{Y}$: Hyperspectral data, $\mathbf{D}$: Endmember library, $\mathbf{r}$: Number of endmembers, $\mu$, $\mu_1$, and $\mu_2$: AL parameters.

**Output:** $\hat{\mathbf{A}}$: Abundances, $\hat{\mathbf{E}}$: Endmembers, $\hat{\mathbf{B}}$: Endmembers' contributions.

**Initialization:** $\mathbf{B} = 0$, $\mathbf{S} = 0$, $\mathbf{L} = 0$, $\mathbf{S}_i = 0$, $\mathbf{L}_i = 0$, $i = 1, 2$

$$
\begin{array}{@{}l@{\quad}l@{}}
\text{A-step :} & \\
\quad \text{for } i = 1 \text{ to } T_1 \text{ do} & \\
\quad \quad \mathbf{A} = \mathrm{QuEC}(\mathbf{A}, \mathbf{S}, \mathbf{L}; \mathbf{Y}, \mathbf{DB}, \mu) \\
\quad \quad \mathbf{S} = \max(0, \mathbf{A} + \mathbf{L}) \\
\quad \quad \mathbf{L} = \mathbf{L} + \mathbf{A} - \mathbf{S} \\
\quad \text{end} & \\
\text{B-step :} & \\
\quad \text{for } i = 1 \text{ to } T_2 \text{ do} & \\
\quad \quad \mathbf{B} = \mathrm{QuEC}(\mathbf{S}_1, \mathbf{L}_1; (\mathbf{S}_2 - \mathbf{L}_2), \mathbf{D}, \mu_1/\mu_2) \\
\quad \quad \mathbf{S}_1 = \max(0, \mathbf{B} + \mathbf{L}_1) \\
\quad \quad \mathbf{S}_2 = (\mathbf{Y}\mathbf{A}^T + \mu_2(\mathbf{D}\mathbf{B} + \mathbf{L}_2))(\mathbf{A}\mathbf{A}^T + \mu_2\mathbf{I}_r)^{-1} \\
\quad \quad \mathbf{L}_1 = \mathbf{L}_1 + \mathbf{B} - \mathbf{S}_1 \\
\quad \quad \mathbf{L}_2 = \mathbf{L}_2 + \mathbf{D}\mathbf{B} - \mathbf{S}_2 \\
\quad \text{end} & \\
\text{end} & \\
\hat{\mathbf{E}} = \mathbf{D}\hat{\mathbf{B}} &
\end{array}
$$

Here, we initialize $\mathbf{S}_1$, $\mathbf{S}_2$, $\mathbf{L}_1$, and $\mathbf{L}_2$ with 0. A-Step and B-Step should be repeated until the convergence otherwise, the cyclic descent with respect to A and B may fail due to the non-convex nature of the problem. The pseudo-code for FaSUn is given in Algorithm 1. Note that for FaSUn the number of endmembers should be given.

### E. SUnS: Sparse Unmixing Using Soft-Shrinkage

Conventional sparse unmixing [15] uses model (3) and sparse regression given by

$$ 
\begin{array}{@{}l@{\quad}l@{}}
\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, \quad \mathbf{1}_m^T \mathbf{X} = \mathbf{1}_n^T,
\end{array}
\quad (22) $$

to estimate the abundances. An ADMM-based algorithm was proposed to solve the problem (22), and therefore, it was called sparse unmixing by variable splitting and augmented Lagrangian (SUnSAL). However, it is suggested to use SUnSAL without ASC due to the conflict with $\ell_1$ [39]. Additionally, ASC was found to be a rigorous constraint that often does not occur in the real world due to noise and signature variability [40]. Therefore, SUnSAL often refers to the problem (22) without ASC. We should note that ignoring ASC breaks physical constraints on pixels for the mixture model. Here, we propose a solution to this challenge. We propose to use an archetypal-type model, i.e., using model (4) and enforce the sparsity on B instead of sum to one. In this way, we can keep ASC while enforcing sparsity. Enforcing $\ell_1$ sparsity on abundance does not allow enforce ASC which is a physical constraint on pixels for the mixture model. ASC does not have to be strict due to noise, however, entirely omitting this constraint may cause abundance estimations that cannot be physically interpreted. Instead, enforcing the sparsity on B allows us to apply ASC on A. Assuming $\mathbf{E} = \mathbf{D}\mathbf{B}$, we enforce the sparsity on matrix B such that only a few atoms of D
contribute to the reconstruction of E. Therefore, we propose
a new optimization,

$$
(\hat{\mathbf{B}}, \hat{\mathbf{A}}) = \arg \min_{\mathbf{B},\mathbf{A}} \frac{1}{2} \|\mathbf{Y} - \mathbf{D}\mathbf{A}\|_{F}^{2} + \lambda \|\mathbf{B}\|_{1} \\
\text{s.t. } \mathbf{A} \ge 0, 1_{r}^{T} \mathbf{A} = 1_{n}^{T}, 0 \le \mathbf{B} \le 1. \tag{23}
$$

One of the main differences of (23) compared to the con-
ventional sparse unmixing is that the former is nonconvex
while the latter is convex if the prior is convex. Moreover,
in the former approach, the number of endmembers must be
predetermined.

Here, we propose an ADMM-based solution to (23). Similar to FaSUn, we use a cyclic descent algorithm, and the **A**-step is the same and therefore we do not repeat it. For the **B**-step, after splitting the variables (**B** = **S**₁ and **DB** = **S**₂) the AL is given by

$$
\arg \min_{\mathbf{B},\mathbf{S}_1,\mathbf{S}_2} \frac{1}{2} \|\mathbf{Y} - \mathbf{S}_2 \mathbf{A}\|_F^2 + \lambda \|\mathbf{S}_1\|_1 + \frac{\mu_1}{2} \|\mathbf{S}_1 - \mathbf{B} - \mathbf{L}_1\|_2^2 \\
+ \frac{\mu_2}{2} \|\mathbf{S}_2 - \mathbf{D}\mathbf{B} - \mathbf{L}_2\|_2^2, \quad \text{s.t.} 0 \le \mathbf{S}_1 \le 1. \tag{24}
$$

The solutions w.r.t. each variable are given by
**B**-Step:

$$
\hat{\mathbf{B}} = (\mu_1 \mathbf{I}_m + \mu_2 \mathbf{D}^T \mathbf{D})^{-1} (\mu_1 (\mathbf{S}_1 - \mathbf{L}_1) + \mu_2 \mathbf{D}^T (\mathbf{S}_2 - \mathbf{L}_2)), \quad (25)
$$

$$
\hat{\mathbf{S}}_1 = \min(1, \max(0, \mathrm{soft}(\mathbf{B} + \mathbf{L}_1, \frac{\lambda}{\mu_1}))), \quad (26)
$$

$$
\hat{\mathbf{S}}_2 = (\mathbf{Y}\mathbf{A}^T + \mu_2(\mathbf{D}\mathbf{B} + \mathbf{L}_2))(\mathbf{A}\mathbf{A}^T + \mu_2\mathbf{I}_r)^{-1}. \quad (27)
$$

$$
\mathbf{L}_1 = \mathbf{L}_1 + \mathbf{B} - \mathbf{S}_1, \qquad (28)
$$

$$
\mathbf{L}_2 = \mathbf{L}_2 + \mathbf{D}\mathbf{B} - \mathbf{S}_2. \qquad (29)
$$

Here, we initialize **S**₁, **S**₂, **L**₁, and **L**₂ with 0. Similarly, **A**-Step and **B**-Step should be repeated until the convergence otherwise, the cyclic descent with respect to **A** and **B** may fail due to the non-convex nature of the problem. The pseudo-code for SUnS is given in Algorithm 2. The initializations are the same as Algorithm 1. Note that for SUnS the number of endmembers should be given.

III. EXPERIMENTAL RESULTS

We employed a total of five datasets, comprising three simu-
lated datasets designed to encompass various mixing scenarios
and two real-world datasets, Houston, and Cuprite, which
is a well-documented geological site. The hyperparameters
for the chosen methods were fine-tuned as per Table I. For
the simulated datasets, we conducted five independent runs,
and the results were then averaged. The standard deviations
are indicated through error bars. We evaluated the perfor-
mance of eight semi-supervised unmixing methods, selected
as follows: SUnSAL [15], CLSUnSAL [17], MUA_SLAVE [20],
S2WSU [19], SUnCNN [22], SUnAA [29], SUnS and Fa-
SUn. The source code used for running SUnSAL, CLSUn-
SAL, MUA_SLAVE, S2WSU, and SUnCNN is available in the
HySUPP toolbox [1] for the sake of reproducibility. Regarding
the parameter tuning procedures, it should be noted that we

**Algorithm 2:** SUnShrink (SUnS)

**Input:** **Y:** Hyperspectral data, **D:** Endmember library, **r:** Number of endmembers, **μ**, **μ**₁, and **μ**₂: AL parameters.

**Output:** **A:** Abundances, **E:** Endmembers, **B:** Endmembers' contributions.

**Initialization:** **B** = 0, **S** = 0, **L** = 0, **S**ᵢ = 0, **L**ᵢ = 0, *i* = 1, 2
Γ = (μ₁**I**<sub>*m*</sub> + μ₂**D**<sup>*T*</sup>**D**)<sup>-1</sup>

for *t* = 1 to *T* do

**A-step:**

for *i* = 1 to *T*<sub>1</sub> do

A = QuEC(A, S, L; Y, DB, μ)

S = max(0, A + L)

L = L + A - S

end

B-step :

for *i* = 1 to *T*<sub>2</sub> do

B = Γ(μ₁(S₁ - L₁) + μ₂D<sup>T</sup>(S₂ - L₂)),

S₁ = min(1, max(0, soft(B + L₁, λ/μ₁))),

S₂ = (YA^T + μ₂(DB + L₂))(AA^T + μ₂I_r)^{-1}

L₁ = L₁ + B - S₁

L₂ = L₂ + DB - S₂

end

end

E = D̂E

tuned the most effective parameters. Indeed, parameter tuning
can be performed for all methods up to some levels. The
conventional sparse unmixing methods were tuned to perform
optimally by tuning their sparsifying regularization parameter.
For SUnCNN, we used the default number of iterations
suggested for different noise levels. SUnAA is a parameter free
technique. FaSUn requires parameter setting for μ₁, μ₂, and
μ₃ since those may affect the convergence. Therefore, those
parameters are set based on the loss function value to ensure
the decrease of the cost function over the iterations. For SUnS
the sparsifying regularization parameter λ was also tuned and
fixed for the simulated and real datasets. Moreover, SUnS
and FaSUn were implemented using PyTorch and provided
in a dedicated Python package called Fast Semisupervised
Unmixing (FUnmix), which is open-source and available at
https://github.com/BehnoodRasti/FUnmix.

In terms of quantitative evaluation, we employed the signal-
reconstruction-error (SRE) measured in decibels (dB) to assess
the estimated abundances, defined by:

$$
\text{SRE}(\mathbf{A}, \hat{\mathbf{A}}) = 20 \log_{10} \frac{\|\mathbf{A}\|_F}{\|\mathbf{A} - \hat{\mathbf{A}}\|_F}. \quad (30)
$$

SRE was chosen over RMSE due to its logarithmic scale
which tends to distinguish the performances better when the
reconstruction error is low which happens often in the case of
simulated datasets. This is a common choice in library-based
unmixing literature [40].

A. Data Description

1) Synthetic *Datasets* with *Spatial Structure*: We simulated two data cubes (DC1 and DC2). DC1 was simulated using a linear mixing model with 5 endmembers selected from the USGS library and 75×75 pixels. The abundance maps are composed of five rows of square regions uniformly distributed over the spatial dimension. This dataset contains pure pixels for all endmembers. DC1 corresponds to the simplest setting where traditional sparse unmixing methods tend to shine. DC2
<table><caption>TABLE I
HYPERPARAMETERS IN DIFFERENT SCENARIOS</caption><thead><tr><th>Methods/Data</th><th>Simulated</th><th>Real</th></tr></thead><tbody><tr><td>SUnSAL</td><td>SNR-dependent (see source code)</td><td>λ = 0.005</td></tr><tr><td>CLSUnSAL</td><td>SNR-dependent (see source code)</td><td>λ = 0.05, μ = 0.01</td></tr><tr><td>MUA_SLAVE</td><td>SNR-dependent (see source code)</td><td>λ<sub>1</sub> = 0.001, λ<sub>2</sub> = 0.01, β = 10, slic_size = 200, slic_reg = 0.01</td></tr><tr><td>S2WSU</td><td>SNR-dependent (see source code)</td><td>λ = 0.001</td></tr><tr><td>SUnCNN</td><td>SNR-dependent (see source code)</td><td>niters = 20000</td></tr><tr><td>SUnAA</td><td>default (see source code)</td><td>default</td></tr><tr><td>SUnS</td><td>T = 10000, T<sub>A</sub> = T<sub>B</sub> = 5, μ<sub>1</sub> = 50, μ<sub>2</sub> = 2, μ<sub>3</sub> = 1, λ = 0.01</td><td>T = 10000, T<sub>A</sub> = T<sub>B</sub> = 5, μ<sub>1</sub> = 400, μ<sub>2</sub> = 100, μ<sub>3</sub> = 1, λ = 0.1</td></tr><tr><td>FaSUn</td><td>T = 10000, T<sub>A</sub> = T<sub>B</sub> = 5, μ<sub>1</sub> = 50, μ<sub>2</sub> = 2, μ<sub>3</sub> = 1</td><td>T = 10000, T<sub>A</sub> = T<sub>B</sub> = 5, μ<sub>1</sub> = 400, μ<sub>2</sub> = 20, μ<sub>3</sub> = 1</td></tr></tbody></table>

has 100×100 pixels and was simulated using a linear mixing model with 9 endmembers. The abundance maps were sampled from a Dirichlet distribution centered at a Gaussian random field to have piece-wise smooth maps with steep transitions. Therefore, DC2 contains spectral variations. DC2 enables us to study the effect of spectral variability in a controlled setting in which the ground truth is available and quantitative results can be obtained. For DC1 and DC2, an endmember library $D \in \mathbb{R}^{224 \times 240}$, composed of 240 spectral signatures was selected from the USGS library with a minimum pair-spectra angle of 4.44°. Synthetic Gaussian noise is added so as to create different signal-to-noise ratio (SNR) scenarios (e.g., 20, 30, and 40 dB SNR).

## 2) Synthetic Datasets with varying Pixel Purity Levels:

We assessed the performance of our chosen methods in an alternative unmixing scenario characterized by the absence of spatial structure but parameterized pixel purity levels. In this context, we were able to explore a spectrum of scenarios ranging from highly mixed, where pure pixels are typically missing, to predominantly pure pixels, with various degrees of mixing in between. The degree of pixel purity was quantified by the parameter $ρ$, with lower values signifying less purity and higher values indicating greater purity.

To construct our dataset, we selected six spectra from the USGS library ($D \in \mathbb{R}^{224 \times 498}$) and created a dataset of size $n = 100 \times 100$ pixels using the following methodology. Initially, we generated a substantial number of abundance samples denoted as $S$. These samples were drawn from the symmetric Dirichlet distribution, employing a scalar concentration parameter $α = 1/r$ where $r$ corresponds to the number of endmembers (in this case, 6). Given a pixel purity level $ρ$, we randomly drew $n$ abundances from $S$ such that their $ℓ_2$ norm fell within the range of $ρ - 0.1$ to $ρ$. Subsequently, we combined the selected spectra based on the sampled abundances to create the final pixel set $Y = [y_1, ..., y_n]$. Synthetic Gaussian noise is eventually added to obtain an input SNR of 30 dB. The objective is to quantify the performance of the proposed methods in a setup where the degree of pixel purity is entirely parametrized. It is used to demonstrate the robustness of our methods in unmixing scenarios that are increasingly harder as $ρ$ decreases.

## 3) Real Dataset with Endmember Bundles:
The Houston dataset used in this paper is a 105×128 pixels subset of an image acquired over the University of Houston campus, Texas, USA, in June 2012. The preprocessed hyperspectral image consists of 144 spectral bands in the 380nm to 1050nm

region [41]. Houston is challenging owing to the spectral variability of the macroscopic materials present in the scene. Therefore the spectral library has been designed in a way to provide several dictionary elements per class. However, there is no library of reference endmembers associated with this image so we use the procedure proposed by [42] to automatically extract endmember bundles from the image pixels. We consider the following four macroscopic materials: concrete, asphalt, vegetation, and red (metallic) roof. In a nutshell, we leverage VCA [2] on ten subsets of the data pixels to estimate four endmembers per subset. We then apply the k-means algorithm (with the cosine similarity measure) using four classes to cluster the estimated endmembers into four endmember bundles. The resulting dictionary consists in 40 endmembers, each associated to a cluster, or bundle.

## 4) Real Dataset with a Spectral Library:
The Cuprite dataset used in this paper contains 250×191 pixels. Cuprite is a well-studied mineral site, and dominant minerals are demonstrated using a geological ground reference. Therefore, the abundance maps estimated by different techniques can be compared visually. We use a library $D \in \mathbb{R}^{188 \times 498}$ composed of 498 spectral pixels from the USGS library. Note that we remove the water absorption and noisy bands, such that the final pixels are of dimension $p = 188$. Cuprite has been selected on account of its challenging characteristics. With almost 500 dictionary elements (from the USGS library), conventional sparse unmixing methods tend to underestimate the abundance of the predominant materials which are Alunite, Chalcedony, and Kaolinite. In addition, due to its high number of pixels (more than 47k), Cuprite offers a testing ground for assessing the efficiency of the proposed algorithms.

## B. Experimental Results: Synthetic Datasets

We first compare the selected methods on the synthetic datasets. Figure 1 summarizes the results in terms of SRE on DC1 and DC2 for different input SNR. The following observations can be formulated:

* In the presence of pure pixels (Fig. 1, DC1), regardless of the noise level, archetypal analysis inspired techniques (i.e. SUnAA, FaSUn, and SUnS) perform very well. MUA_SLAVE is suitable when the noise is significant (e.g. 20 dB) but its performance drops when dealing with less noisy images due to its segmentation approach. CLSUnSAL performs a bit better than SUnSAL but struggles to compete with the top performing methods overall.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Abundance SRE (↑) in dB for the selected semi-supervised methods on two simulated datasets (DC1 and DC2) for three noise levels.</figcaption></figure>

S2WSU and SUnCNN provide similar results on this dataset. It is worth mentioning that most methods, except the archetypal analysis inspired ones, require tuning their regularization parameter to obtain competitive results depending on the input SNR, which is a major hindrance.

* In the presence of spectral variability (Fig. 1, DC2), SUnCNN performs very well, likely due to its convolutional architecture that is suited to capture the spectral variability relying on the spatial structure of the data. Similarly, S2WSU obtains competitive results. SUnAA and FaSUn, which both solve problem (5), outperform the other methods. It should be noted that SUnS struggles

in the lower SNR scenarios (i.e. 20 and 30 dB) which reveals the advantage of convexity constraint compared to the sparsity prior. Moreover, it appears that SUnSAL and CLSUnSAL demonstrate the poorest performances compared to the other methods.

Figure 2 summarizes the results in terms of SRE for different pixel purity levels using a fixed SNR (30 dB). The following observations can be made:

* The archetypal analysis inspired models severely outperform their sparse unmixing counterparts, regardless of the pixel purity level. This is particularly striking when the pixel purity is low, meaning the image only
contains highly mixed pixels. This indicates that SUnAA, FaSUn and SUnS are better suited to handle highly mixed scenarios. Furthermore, the latter methods do not leverage spatial information contrarily to S2WSU, whose performance drops significantly due to the absence of spatial structure in the data.

* It is worth mentioning that the endmembers library, D, has not been pruned, contrarily to the previous setups, in which the number of atoms in the dictionary went down from 498 to 240. Therefore there is a clear benefit in having access to the number of endmembers present in the scene (i.e. *r*), which is available for SUnAA, FaSUn and SUnShrink.

Figure 3 and 4 demonstrate the visual comparisons of estimated abundances by applying different semisupervised unmixing techniques to DC1 and DC2, respectively, for SNR=20 dB. Overall, the visual comparisons confirm that FaSUn and SUnA perform similarly and outperform the other semisupervised techniques for those datasets. SUnS performs similarly to SUnAA and FaSun in the case of DC1. However, in the case of DC2, SUnS cannot successfully estimate the abundance map associated with endmember 7. It is worth mentioning that, in the case of DC1 (20 dB), MUA_SLAVE provides the highest SRE, however, the visual comparisons reveal that abundances are oversmoothed due to the prior segmentation step which can be associated with the high SRE only for low SNR. SUnCNN performs well in the case of DC2 but the abundances for DC1 are oversmoothened. The abundances estimated by SUnSAL, CLSUnSAL, and S2WSU are not competitive with the other methods.

### C. Experimental Results: Real Data

Figure 5 compares the estimated abundances by applying different semisupervised unmixing methods to the Houston dataset. In this experiment, we use endmember bundles as a library. The abundances are shown by applying a sum to one of the associated abundances to endmember bundles except for SUnAA, SUnS, and FaSUn since they use ASC. The results of SUnAA and FaSUn are similar since they solve the same problem. FaSUn, SUnAA, and SUnCNN demonstrate a higher abundance ratio for Concrete compared to the other techniques. SUnS and SUnSAL provide higher abundance ratios for Asphalt and Red Roof, respectively. The sparse unmixing techniques using the SR model provide a higher abundance ratio for vegetation. Overall, FaSUn SUnAA, and SUnCNN show better visual performance in distinguishing the materials.

Figure 6 visually compares the estimated abundances for three dominant materials, i.e., Chalcedony, Alunite and Kaolinite, using the geological map as a reference. It is worth mentioning that hyperparameters for each method had to be tuned, except for SUnAA as it is parameter-free. Moreover, for the archetypal analysis-inspired methods to work, the number of endmembers in the scene (i.e. *r*) was set to *r* = 14.

Visual comparison based on the reference map reveals that SUnAA better estimated Chalcedony compared to the other methods. FaSUn does not exhibit the same saliency as SUnAA,

but still detects Chalcedony on a bigger area of the map than other sparse unmixing methods. As for Alunite, SUnAA and FaSUn both show strong responses to the mineral in the expected areas. Sharp abundance maps are obtained for Kaolinite by SUnAA and FaSUn which are in line with the reference map. Overall, SUnS shows similar performances as the other sparse unmixing techniques. It is worth mentioning that SUnS introduces another regularization parameter, λ, (similar to the other sparse unmixing methods) that requires additional tuning, compared to FaSUn.

### D. Processing Time and Computational Complexity

Perhaps the gist of our contributions is the considerable scalability of our proposed approaches (FaSUn and SUnShrink), as highlighted in Table II. Note that SUnCNN processing time depends on the number of iterations, which is itself dependent on the input SNR. Here we report the processing time for a fixed SNR equal to 30 dB obtained using a computer with an Intel(R) Xeon(R) Silver 4110 CPU at 2.10GHz, 32 cores, 64 Gb of RAM, and a NVIDIA GeForce RTX 2080 Ti GPU with 12 Gb of RAM. We observe that SUnAA processing time is prohibitive for real-world applications, as soon as the number of pixels exceeds *n* = 10000. On the contrary, the most efficient methods are SUnS and FaSUn on Cuprite, which exploit GPU computations despite using many outer iterations, i.e. *T* = 10000. It is worth highlighting the growth rate of consumed time by the FaSUn and SUnS compared to the other techniques. They take around three minutes to perform on a dataset with 90k pixels and 224 bands. These results shed light on the scalability challenge posed to semi-supervised methods despite ongoing efforts. We should note that among the methods mentioned above, SUnCNN, FaSUn, and SUnS use GPU and the rest use CPU for computations. Concerning the computational complexity, each iteration of FaSUn and SUnS algorithms is $O((T_1 + T_2)npr)$ due to costly matrix products and solving linear systems.

## IV. CONCLUSION

We proposed two nonconvex optimizations for hyperspectral unmixing relying on a new type of linear model called FaSUn and SUnS. FaSun utilizes convexity constraint while SUnS uses a sparse prior. We derived ADMM-based solutions for those problems implemented using PyTorch. We conducted a comprehensive evaluation of our proposed techniques by comparing them with SOTA methods. This evaluation was performed on three simulated datasets, considering spatial structure, spectral variability, and various pixel purity and noise levels, and we measured their performance in terms of signal-to-reconstruction error. The results strongly support the superiority of the proposed unmixing technique, FaSUn, over existing state-of-the-art methods. Notably, FaSUn consistently achieved better performance across a range of scenarios. The results confirmed the advantage of the convexity constraint compared to the sparsity-promoting prior. Additionally, we applied these unmixing techniques to real-world data using the Cuprite dataset. To validate the accuracy and practical utility of the proposed methods, we visually compared the results
<figure><img src="image_3.png" /><figcaption>Fig. 2. Abundance SRE (↑) in dB for the selected semi-supervised methods on three different pixel purity levels (ρ) for a given input SNR (30 dB)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Visual comparisons of abundance maps estimated by using different semi-supervised unmixing methods applied to DC1 (20 dB).</figcaption></figure>

with geological reference maps. Furthermore, the proposed ADMM-based algorithms demonstrated remarkable efficiency in addition to their superior performance. This was evident when comparing the processing times for large datasets, highlighting the practical advantages of our approach.

## APPENDIX A
### DERIVATION OF QUEC FUNCTION

Assuming the quadratic programming (or least squares) with the equality constraint as

$$ \begin{gathered} \hat{\mathbf{A}} = \arg \min_{\mathbf{A}} \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \frac{\mu}{2} \|\mathbf{S} - \mathbf{A} - \mathbf{L}\|_F^2 \\ \text{s.t.} \quad \mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T, \end{gathered} \quad (31) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Visual comparisons of abundance maps estimated by using different semi-supervised unmixing methods applied to DC2 (20 dB).</figcaption></figure>

the Lagrangian function is given by

$$ \mathcal{L}(\mathbf{A}, \nu) = \frac{1}{2} \|\mathbf{Y} - \mathbf{E}\mathbf{A}\|_F^2 + \frac{\mu}{2} \|\mathbf{S} - \mathbf{A} - \mathbf{L}\|_F^2 + \nu^T (\mathbf{1}_r^T \mathbf{A} - \mathbf{1}_n^T), \quad (32) $$

where the solution can be given using Karush-Kuhn-Tucker (KKT) conditions. Here, we derive the KKT conditions for (32). The stationary condition is given by $\nabla_{\mathbf{A}}\mathcal{L}(\mathbf{A}, \nu) = 0$,

$$ (\mathbf{E}^T \mathbf{E} + \mu \mathbf{I}) \mathbf{A} + \mathbf{1}_r \nu = \mathbf{E}^T \mathbf{Y} + \mu (\mathbf{S} - \mathbf{L}) \quad (33) $$

For the primal feasibility we hold $\mathbf{1}_r^T \mathbf{A} = \mathbf{1}_n^T$. Therefore, we have

$$ \begin{pmatrix} \mathbf{E}^T \mathbf{E} + \mu \mathbf{I} & \mathbf{1}_r \\ \mathbf{1}_r^T & 0 \end{pmatrix} \begin{pmatrix} \mathbf{A} \\ \nu \end{pmatrix} = \begin{pmatrix} \mathbf{E}^T \mathbf{Y} + \mu (\mathbf{S} - \mathbf{L}) \\ \mathbf{1}_n^T \end{pmatrix}, \quad (34) $$

using the blockwise inversion, the solution is given by

$$ \hat{\mathbf{A}} = (\mathbf{Q} + \mathbf{Q}\mathbf{1}_r c\mathbf{1}_r^T \mathbf{Q})(\mathbf{E}^T \mathbf{Y} + \mu(\mathbf{S} - \mathbf{L})) - \mathbf{Q}\mathbf{1}_r c\mathbf{1}_n^T \quad (35) $$

Where

$$ \mathbf{Q} = (\mathbf{E}^T \mathbf{E} + \mu \mathbf{I}_r)^{-1} \quad (36) $$

$$ c = -1/(\mathbf{1}_r^T \mathbf{Q} \mathbf{1}_r). \quad (37) $$

As can be seen, the augmented term in (32), makes matrix Q to be always non-singular (note that $\mu > 0$) and therefore the closed form solution (35) becomes feasible.

## ACKNOWLEDGMENT

AZ, JM, and JC were supported by ANR 3IA MIAI@Grenoble Alpes (ANR-19-P3IA-0003). AZ and JM were supported by the ERC grant number 101087696 (APHELAIA project). We would like to thank the Hyperspectral Image Analysis Laboratory at the University of Houston, and the IEEE GRSS Image Analysis and Data
<figure><img src="image_6.png" /><figcaption>Fig. 5. Estimated abundances obtained by different semi-supervised methods on the Houston dataset.</figcaption></figure>

<table><caption>TABLE II
PROCESSING TIME. THE BEST RESULTS ARE IN BOLDFACE AND THE SECOND BEST ARE UNDERLINED.</caption><thead><tr><th></th><th># Pixels</th><th># Bands</th><th># Endmembers</th><th># Atoms</th><th>SUnSAL</th><th>CLSUnSAL</th><th>MUA_SLAVE</th><th>S2WSU</th><th>SUnCNN</th><th>SUnAA</th><th>SUnS</th><th>FaSUn</th></tr></thead><tbody><tr><td>DC1</td><td>5625</td><td>224</td><td>5</td><td>240</td><td><strong>7.1</strong></td><td>15.3</td><td><strong>6.1</strong></td><td>41.1</td><td>94.7</td><td>131.1</td><td>85.3</td><td>89.2</td></tr><tr><td>DC2</td><td>10000</td><td>224</td><td>9</td><td>240</td><td><strong>18.9</strong></td><td>27.2</td><td><strong>10.2</strong></td><td>51.0</td><td>78.0</td><td>264.8</td><td>80.0</td><td>77.6</td></tr><tr><td>Mixed Pixels 10k</td><td>10000</td><td>224</td><td>6</td><td>498</td><td><strong>46.4</strong></td><td>159.1</td><td>98.2</td><td><strong>61.3</strong></td><td>87.2</td><td>208.7</td><td>136.9</td><td>134.3</td></tr><tr><td>Houston</td><td>13440</td><td>144</td><td>4</td><td>40</td><td><strong>10.9</strong></td><td>46.4</td><td><strong>5.1</strong></td><td>28.4</td><td>220.4</td><td>109.6</td><td>61.8</td><td>63.1</td></tr><tr><td>Cuprite</td><td>47750</td><td>188</td><td>14</td><td>498</td><td>200.3</td><td>343.8</td><td>198.5</td><td>662.0</td><td>769.6</td><td>1838.0</td><td><strong>154.5</strong></td><td><strong>147.7</strong></td></tr><tr><td>Mixed Pixels 90k</td><td>90000</td><td>224</td><td>6</td><td>498</td><td>566.5</td><td>2065.3</td><td>755.9</td><td>660.0</td><td>576.0</td><td>1359.0</td><td><strong>183.4</strong></td><td><strong>194.6</strong></td></tr></tbody></table>

Fusion Technical Committee for providing the Houston dataset.

REFERENCES

[1] B. Rasti, A. Zouaoui, J. Mairal, and J. Chanussot, "Image processing and machine learning for hyperspectral unmixing: An overview and the hysupp python package," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 62, pp. 1–31, 2024.

[2] J. Nascimento and J. Bioucas-Dias, "Vertex component analysis: A fast algorithm to extract endmembers spectra from hyperspectral data," in *Pattern Recognition and Image Analysis*, Francisco José Perales, Aurélio J. C. Campilho, Nicolás Pérez de la Blanca, and Alberto Sanfeliu, Eds., Berlin, Heidelberg, 2003, pp. 626–635, Springer Berlin Heidelberg.

[3] R. Heylen, D. Burazerovic, and P. Scheunders, "Fully constrained least squares spectral unmixing by simplex projection," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 49, no. 11, pp. 4112–4122, Nov 2011.

[4] J. Li and J. M. Bioucas-Dias, "Minimum volume simplex analysis: A fast algorithm to unmix hyperspectral data," in *IGARSS 2008 - 2008 IEEE International Geoscience and Remote Sensing Symposium*, 2008, vol. 3, pp. III – 250–III – 253.

[5] J. M. Bioucas-Dias, "A variable splitting augmented lagrangian approach to linear spectral unmixing," in *2009 First Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing*, 2009, pp. 1–4.

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

[7] C.-I. Chang, X.-L. Zhao, M.L.G. Althouse, and J. J. Pan, "Least squares subspace projection approach to mixed pixel classification for hyperspectral images," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 36, no. 3, pp. 898–912, 1998.

[8] J.J. Settle and N. Drake, "Linear mixing and the estimation of ground cover proportions," *International Journal of Remote Sensing - INT J REMOTE SENS*, vol. 14, pp. 1159–1177, 04 1993.

[9] C.-I. Chang and D. C. Heinz, "Constrained subpixel target detection for remotely sensed imagery," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 38, no. 3, pp. 1144–1159, 2000.

[10] Y.E. Shimabukuro, *Shade Images Derived from Linear Mixing Models of Multispectral Measurements of Forested Areas*, Colorado State University, 1987.

[11] Y.E. Shimabukuro and J.A. Smith, "The least-squares mixing models to generate fraction images derived from remote sensing multispectral data," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 29, no. 1, pp. 16–20, 1991.

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

[13] J. Eckstein and D. P. Bertsekas, "On the douglas-rachford splitting
<figure><img src="image_7.png" /><figcaption>(a) Reference Map</figcaption></figure>

<caption>(b) Estimated abundances</caption>

<caption>Fig. 6. Estimated abundances obtained by applying different semi-supervised methods to Cuprite compared with the geological reference map.</caption>

method and the proximal point algorithm for maximal monotone operators,” Mathematical Programming, vol. 55, pp. 293–318, 1992.

[14] C. Quintano, A. Fernandez-Manso, and D. A. Roberts, “Multiple endmember spectral mixture analysis (mesma) to map burn severity levels from landsat images in mediterranean countries,” *Remote Sensing of Environment*, vol. 136, pp. 76–88, 2013.

[15] J. M. Bioucas-Dias and M. A. T. Figueiredo, “Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing,” in *2nd Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing*, 2010, pp. 1–4.

[16] M. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Total variation spatial regularization for sparse hyperspectral unmixing,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 50, no. 11, pp. 4484–4502, 2012.

[17] M. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Collaborative sparse regression for hyperspectral unmixing,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 52, no. 1, pp. 341–354, 2014.

[18] R. Wang, H.-C. Li, A. Pizurica, J. Li, A. Plaza, and W. J. Emery, “Hyperspectral unmixing using double reweighted sparse regression and total variation,” *IEEE Geoscience and Remote Sensing Letters*, vol. 14, no. 7, pp. 1146–1150, 2017.

[19] S. Zhang, J. Li, H. Li, C. Deng, and A. Plaza, “Spectral–spatial weighted sparse regression for hyperspectral image unmixing,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 56, no. 6, pp. 3265–3276, 2018.

[20] R. A. Borsoi, T. Imbiriba, J. C. M. Bermudez, and C. Richard, “A fast multiscale spatial regularization for sparse hyperspectral unmixing,” *IEEE Geoscience and Remote Sensing Letters*, vol. 16, no. 4, pp. 598–602, 2019.

[21] T. Ince, “Superpixel-based graph laplacian regularization for sparse hyperspectral unmixing,” *IEEE Geoscience and Remote Sensing Letters*, pp. 1–5, 2020.

[22] B. Rasti and B. Koirala, “SUnCNN: Sparse unmixing using unsupervised convolutional neural network,” *IEEE Geoscience and Remote Sensing Letters*, vol. 19, pp. 1–5, 2022.

[23] Y. Lin and P. Gader, “Addressing spectral variability in hyperspectral unmixing with unsupervised neural networks,” in *2022 12th Workshop on Hyperspectral Imaging and Signal Processing: Evolution in Remote Sensing (WHISPERS)*, 2022, pp. 1–5.

[24] Q. Qian, F. Xiong, and J. Zhou, “Deep unfolded iterative shrinkage-thresholding model for hyperspectral unmixing,” in *IGARSS 2019* -

2019 IEEE International Geoscience and Remote Sensing Symposium, 2019, pp. 2151–2154.

[25] Y. Shao, Q. Liu, and L. Xiao, “Iviu-net: Implicit variable iterative unrolling network for hyperspectral sparse unmixing,” *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, vol. 16, pp. 1756–1770, 2023.

[26] Y. Qian, F. Xiong, Q. Qian, and J. Zhou, “Spectral mixture model inspired network architectures for hyperspectral unmixing,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 58, no. 10, pp. 7418–7434, 2020.

[27] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” *Communications on Pure and Applied Mathematics*, vol. 57, no. 11, pp. 1413–1457, 2004.

[28] C. Zhou and M. R. D. Rodrigues, “Admm-based hyperspectral unmixing networks for abundance and endmember estimation,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 60, pp. 1–18, 2022.

[29] B. Rasti, A. Zouaoui, J. Mairal, and J. Chanussot, “Sunaa: Sparse unmixing using archetypal analysis,” *IEEE Geoscience and Remote Sensing Letters*, pp. 1–1, 2023.

[30] Y. Chen, J. Mairal, and Z. Harchaoui, “Fast and robust archetypal analysis for representation learning,” *2014 IEEE Conference on Computer Vision and Pattern Recognition*, pp. 1478–1485, 2014.

[31] J.M. Bioucas-Dias and J.M.P. Nascimento, “Hyperspectral subspace identification,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 46, no. 8, pp. 2435–2445, 2008.

[32] B. Rasti, M.O. Ulfarsson, and J.R. Sveinsson, “Hyperspectral subspace identification using SURE,” *IEEE Geoscience and Remote Sensing Letter*, vol. 12, no. 12, pp. 2481–2485, Dec. 2015.

[33] C.-I Chang and Q. Du, “Estimation of number of spectrally distinct signal sources in hyperspectral imagery,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 42, no. 3, pp. 608–619, 2004.

[34] J. Harsanyi, W. Farrand, and C.I Chang, “Determining the number and identity of spectral endmembers: An integrated approach using neyman?pearson eigenthresholding and iterative constrained rms error minimization,” in *9th Thematic Conf. Geologic Remote Sens.*, 1993.

[35] B. Luo, J. Chanussot, S. Doute, and L. Zhang, “Empirical automatic estimation of the number of endmembers in hyperspectral images,” *IEEE Geoscience and Remote Sensing Letters*, vol. 10, no. 1, pp. 24–28, 2013.

[36] M. Berman, H. Kiiveri, R. Lagerstrom, A. Ernst, R. Dunne, and J.F. Huntington, “Ice: a statistical approach to identifying endmembers in
hyperspectral images,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 42, no. 10, pp. 2085–2095, 2004.

[37] A. Zare and P. Gader, “Sparsity promoting iterated constrained endmember detection in hyperspectral imagery,” *IEEE Geoscience and Remote Sensing Letters*, vol. 4, no. 3, pp. 446–450, 2007.

[38] A Ambikapathi, T.-H. Chan, C.-Y. Chi, and K. Keizer, “Hyperspectral data geometry-based estimation of number of endmembers using p-norm-based pure pixel identification algorithm,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 51, no. 5, pp. 2753–2769, May 2013.

[39] J. Sigurdsson, M. O. Ulfarsson, and J. R. Sveinsson, “Hyperspectral unmixing with $l_q$ regularization,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 52, no. 11, pp. 6793–6806, 2014.

[40] M. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Sparse unmixing of hyperspectral data,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 49, no. 6, pp. 2014–2039, 2011.

[41] C. Debes, A. Merentitis, R. Heremans, J. Hahn, N. Frangiadakis, T. van K., W. Liao, R. Bellens, A. Pižurica, S. Gautama, W. Philips, S. Prasad, Q. Du, and F. Pacifici, “Hyperspectral and lidar data fusion: Outcome of the 2013 grss data fusion contest,” *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, vol. 7, no. 6, pp. 2405–2418, 2014.

[42] B. Somers, M. Zortea, A. Plaza, and G. P. Asner, “Automated extraction of image-based endmember bundles for improved spectral unmixing,” *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, vol. 5, no. 2, pp. 396–408, 2012.

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

**Julien Mairal** received the graduate degree from Ecole Polytechnique, Palaiseau, France, in 2005, and the Ph.D. degree from Ecole Normale Supérieure, Cachan, France, in 2010. After that, he joined the statistics department at UC Berkeley as a post-doctoral researcher. In 2012, he joined Inria, Grenoble, France, where he is currently a research director and head of the Thoth team. His research interests include machine learning, computer vision, mathematical optimization, and statistical image and signal processing. He received a starting grant and a

consolidator grant from the European Research Council, respectively in 2016 and 2022. He was awarded the Cor Baayen prize in 2013, the IEEE PAMI young research award in 2017 and the test-of-time award at ICML 2019.

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

**Behnood Rasti (M'12-SM'19)** received the B.Sc. and M.Sc. degrees in electronics and electrical engineering from the Electrical Engineering Department, University of Guilan, Rasht, Iran, in 2006 and 2009, respectively, and the Ph.D. degree in electrical and computer engineering from the University of Iceland, Reykjavik, Iceland, in 2014. He was a Valedictorian as an M.Sc. Student in 2009. In 2015 and 2016, he worked as a Post-Doctoral Researcher with the Electrical and Computer Engineering Department, University of Iceland. From 2016 to 2019,

he was a Lecturer with the Center of Engineering Technology and Applied Sciences, Department of Electrical and Computer Engineering, University of Iceland. Dr. Rasti was a Humboldt Research Fellow in 2020 and 2021. From 2022 to 2023, he was a Principal Research Associate with Helmholtz-Zentrum Dresden-Rossendorf (HZDR), Dresden, Germany. He is currently a senior research scientist at the Faculty of Electrical Engineering and Computer Science, Technische Universität Berlin, and the Berlin Institute for the Foundations of Learning and Data, Berlin, Germany. His research interests include signal and image processing, machine/deep learning, remote sensing, and artificial intelligence.

Dr. Rasti won the Doctoral Grant of the University of Iceland Research Fund “The Eimskip University Fund” and the “Alexander von Humboldt Research Fellowship Grant” in 2013 and 2019, respectively. He serves as an Associate Editor for the IEEE GEOSCIENCE AND REMOTE SENSING LETTERS (GRSL).

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

**Jocelyn Chanussot (M'04-SM'04-F'12)** received the M.Sc. degree in electrical engineering from the Grenoble Institute of Technology (Grenoble INP), Grenoble, France, in 1995, and the Ph.D. degree from the Université de Savoie, Annecy, France, in 1998. From 1999 to 2023, he has been with Grenoble INP, where he was a Professor of signal and image processing. He is currently a Research Director with INRIA, Grenoble. His research interests include image analysis, hyperspectral remote sensing, data fusion, machine learning and artificial intelligence.

He has been a visiting scholar at Stanford University (USA), KTH (Sweden) and NUS (Singapore). Since 2013, he is an Adjunct Professor of the University of Iceland. In 2015-2017, he was a visiting professor at the University of California, Los Angeles (UCLA). He holds the AXA chair in remote sensing and is an Adjunct professor at the Chinese Academy of Sciences, Aerospace Information research Institute, Beijing. Dr. Chanussot is the founding President of IEEE Geoscience and Remote Sensing French chapter (2007-2010) which received the 2010 IEEE GRS-S Chapter Excellence Award. He has received multiple outstanding paper awards. He was the Vice-President of the IEEE Geoscience and Remote Sensing Society, in charge of meetings and symposia (2017-2019). He was the General Chair of the first IEEE GRSS Workshop on Hyperspectral Image and Signal Processing, Evolution in Remote sensing (WHISPERS). He was the Chair (2009-2011) and Cochair of the GRS Data Fusion Technical Committee (2005-2008). He was a member of the Machine Learning for Signal Processing Technical Committee of the IEEE Signal Processing Society (2006-2008) and the Program Chair of the IEEE International Workshop on Machine Learning for Signal Processing (2009). He is an Associate Editor for the IEEE Transactions on Geoscience and Remote Sensing, the IEEE Transactions on Image Processing and the Proceedings of the IEEE. He was the Editor-in-Chief of the IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing (2011-2015). In 2014 he served as a Guest Editor for the IEEE Signal Processing Magazine. He is a Fellow of the IEEE, an ELLIS Fellow, a Fellow of the Asia-Pacific Artificial Intelligence Association, a member of the Institut Universitaire de France (2012-2017) and a Highly Cited Researcher (Clarivate Analytics/Thomson Reuters, since 2018).

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

**Alexandre Zouaoui** received his M.Sc. degree in computer science from Télécom Paris, Paris, France, in 2019. He received his Ph.D. in applied mathematics at Université Grenoble Alpes, Grenoble, France, in 2024 on Sparse and Archetypal Decomposition Algorithms for Hyperspectral Images Restoration and Spectral Unmixing. His research interests include computer vision, interpretable machine learning and remote sensing image processing. In 2024, he joined Data Science Experts as a research engineer focusing on applications related to hyperspectral imagery.