Learning Interpretable Deep Disentangled Neural Networks for
  Hyperspectral Unmixing
Although considerable effort has been dedicated to improving the solution to
the hyperspectral unmixing problem, non-idealities such as complex radiation
scattering and endmember variability negatively impact the performance of most
existing algorithms and can be very challenging to address. Recently, deep
learning-based frameworks have been explored for hyperspectral umixing due to
their flexibility and powerful representation capabilities. However, such
techniques either do not address the non-idealities of the unmixing problem, or
rely on black-box models which are not interpretable. In this paper, we propose
a new interpretable deep learning method for hyperspectral unmixing that
accounts for nonlinearity and endmember variability. The proposed method
leverages a probabilistic variational deep-learning framework, where
disentanglement learning is employed to properly separate the abundances and
endmembers. The model is learned end-to-end using stochastic backpropagation,
and trained using a self-supervised strategy which leverages benefits from
semi-supervised learning techniques. Furthermore, the model is carefully
designed to provide a high degree of interpretability. This includes modeling
the abundances as a Dirichlet distribution, the endmembers using
low-dimensional deep latent variable representations, and using two-stream
neural networks composed of additive piecewise-linear/nonlinear components.
Experimental results on synthetic and real datasets illustrate the performance
of the proposed method compared to state-of-the-art algorithms.
Ricardo Augusto Borsoi, Deniz Erdoğmuş, Tales Imbiriba
# Learning Interpretable Deep Disentangled Neural Networks for Hyperspectral Unmixing

###### Abstract

Although considerable effort has been dedicated to improving the solution to the hyperspectral unmixing problem, non-idealities such as complex radiation scattering and endmember variability negatively impact the performance of most existing algorithms and can be very challenging to address. Recently, deep learning-based frameworks have been explored for hyperspectral unmixing due to their flexibility and powerful representation capabilities. However, such techniques either do not address the non-idealities of the unmixing problem, or rely on black-box models which are not interpretable. In this paper, we propose a new interpretable deep learning method for hyperspectral unmixing that accounts for nonlinearity and endmember variability. The proposed method leverages a probabilistic variational deep-learning framework, where disentanglement learning is employed to properly separate the abundances and endmembers. The model is learned end-to-end using stochastic backpropagation, and trained using a self-supervised strategy which leverages benefits from semi-supervised learning techniques. Furthermore, the model is carefully designed to provide a high degree of interpretability. This includes modeling the abundances as a Dirichlet distribution, the endmembers using low-dimensional deep latent variable representations, and using two-stream neural networks composed of additive piecewise-linear/nonlinear components. Experimental results on synthetic and real datasets illustrate the performance of the proposed method compared to state-of-the-art algorithms.

 Hyperspectral data, spectral unmixing, neural networks, disentanglement, deep learning.

## I Introduction

### _Background_

Due to physical limitations of the imaging process, hyperspectral images (HIs) provide very high spectral resolution but low spatial resolution, which means that a pixel usually contains a mixture of several different materials [2]. Hyperspectral unmixing (HU) consists in estimating the spectral signatures of pure materials (i.e., _endmembers_ - EMs) in a scene and the proportions with which they are contained in each pixel (i.e., _abundances_) directly from an HI [3]. Due to the unsupervised nature of HU, adequately exploring the physics of the problem when devising modeling strategies is paramount for obtaining stable and high-quality EM and abundance estimations. Simplistic methods considered the interaction between light and the EMs to be linear [3]. However, this model is over-simplified, degrading the quality of the estimates. Thus, addressing important non-idealities such as nonlinear interactions between light and the materials [4] and the variability of the EMs in different HI pixels [5] has become the subject of much attention more recently.

Different approaches have been proposed to perform HU (see Section II). However, traditional methods often lack the flexibility to represent non-idealities observed in practical HIs. This motivated the use of machine learning approaches for HU that ally both flexibility and performance [6, 7]. Nonetheless, interpretability remains a key point when leveraging machine learning strategies in HI analysis [8].

Recently, physically motivated machine learning approaches have been successfully applied to HU [9, 10, 11, 12]. The advantage of such models with respect to fully black-box strategies lies in the interpretability of the estimated EMs and abundance parameters, which is a requirement for meaningful unmixing results as it allows a user to understand the causes for the algorithm's behavior [13]. When deep learning strategies come into play, autoencoder (AEC) architectures are of special interest due to the intrinsic low-dimensionality of the abundance space with respect to the pixels, and to the connection between such strategies and hyperspectral mixing models [14, 15]. Thus, several approaches using AECs were proposed to solve HU addressing phenomena such as nonlinearity [10, 16, 17], EM variability [11, 18] and outliers [19].

Although such deep learning methods presented relevant solutions for HU reaching high levels of accuracy while retaining physical interpretation, such strategies fail to provide a separation between EMs and abundances that is both interpretable and accounts for existing spectral variability and nonlinear effects. Recently, supervised disentanglement learning has become a popular approach to separate latent variables in deep learning models into different factors of variation that can have a physical interpretation [20]. Disentangled decompositions have been considered for different applications (e.g., separating content from style in images [20]), and its potential will be explored in this work to aid the separation between abundance and EM variations in HU.

### _Contributions of this work_

In this work, we aim to develop a learning-based unmixing algorithm which accounts for challenges including both nonlinearities in the mixing process and spectral variability of the EMs. Moreover, different from black-box models, the
###### Abstract

In this paper, we propose a novel framework for multi-dimensional Gaussian random variable (EM), which combines the parameters of the model and an approximate posterior distribution. This approach is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable. We propose a novel framework for multi-dimensional Gaussian random variable (EM), which combines the parameters of the model and an approximate posterior distribution. This approach is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM). We propose a novel framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional random variable (EM), which is designed to provide a framework for multi-dimensional random variable (EM), which is designed to provide a framework for multi-dimensional random variable (EM), which is designed to provide a framework for multi-dimensional random variable (EM), which is designed to provide a framework for multi-dimensional random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional random variable (EM), which is designed to provide a framework for multi-dimensional Gaussian random variable (EM), which is designed to provide a framework for multi-dimensional random variable (EM), which is designed to provide
[52], combinations of additive and multiplicative factors [53], or multiscale spatial models [54]. HU is then formulated as a non-convex optimization problem using carefully designed regularizations to constrain the degrees of freedom of the model. Recent work also considered a linear-quadratic non-linear mixture model with EM variability in HU [55].

Deep learning-based HUDeep learning has become an established tool for hyperspectral imaging, with methods such as graph convolutional neural networks (NNs) [56] and the transformer architecture [57] bringing significant advances to hyperspectral classification. Various deep learning-based approaches have also been proposed for HU. Early works trained NNs as supervised regression methods to learn a mapping from the mixed pixels to their abundance values [58, 59]. However, those methods are hindered by the limited availability of training data with ground truth. This motivated the development of algorithms such as extended support vector machines (i.e., soft classification) [60], which require only EM signatures as training samples.

Recently, autoencoder (AEC) networks have become widely used in unsupervised HU [7]. AECs are encoder-decoder (i.e., bottleneck) NNs which can be applied for HU by identifying the decoder with the mixing model (1), the encoder with its inverse, and the latent codes with the abundances [14, 15]. AECs have been developed for linear HU using, e.g., appropriately designed encoder networks or preprocessing approaches to reduce the effect of noise and outliers [61, 19], sparsity constraints [62, 63], or convolutional architectures to exploit spatial information about the HIs [64, 65]. Coupled convolutional AECs were also used to jointly unmix hyperspectral and multispectral data for image fusion [66].

Other methods based on AECs proposed to address non-linear HU by using nonlinear architectures for the decoder, including a post-nonlinear model [16], additive nonlinearities [67, 68], or using specifically designed nonlinear NN layers [69]. A model-based architecture for the encoder was also proposed in [17] by exploring the relationship between the AEC and the nonlinear mixing model.

AECs have also been used to address EM variability in HU as generative EM models to capture the low-dimensional manifold of spectral variability [11]. This model was used in matrix factorization [11], structured sparse regression [21], and probabilistic [18] HU methods. Endmember variability was also addressed in [70] a spatial-spectral autoencoder using an additive perturbation model [50], and in non-AEC-based approaches using Gaussian Process regression [71, 72]. Sets of multiple pure pixels extracted from the HI (which can be seen as samples from EM signatures) have also been explored to improve the robustness of AECs in HU, either by learning mappings from the pure pixels to their corresponding abundances [10, 73], or by regularizing the EMs estimated by the decoder [74]. Recurrent NNs have also been recently employed in [75] to perform HU accounting for the temporal variability of the EMs.

Other works have also considered different loss functions for AEC training, such as the Wasserstein distance [76] and adversarial losses [77]. A cycle-consistency loss was used in [78] to guide the reconstruction of two cascaded AECs in HU. Self-supervised learning has also been integrated into the design of AEC strategies for HU [79]. Moreover, non-AEC approaches were also proposed for HU, such as parametrizations of the abundances using untrained deep prior models [80] or designing NNs for unmixing by unfolding optimization algorithms such as the ADMM [81] or sparse regression methods [12]. Differently from approaches such as unfolded NNs, the proposed method is built upon a statistical model of the mixing and unmixing process, where some of the PDFs are parametrized by learnable models.

## III Proposed method

DefinitionsLet us consider a dataset with \(N_{U}\) unlabeled HI pixels \(\mathcal{D}_{U}=\{\mathbf{y}_{1},\ldots,\mathbf{y}_{N_{U}}\}\). We also consider a supervised dataset \(\mathcal{D}_{S}=\{(\mathbf{y}_{1},\mathbf{a}_{1},\mathbf{M}_{1}),\ldots,(\mathbf{y}_{N_{S}}, \mathbf{a}_{N_{S}},\mathbf{M}_{N_{S}})\}\) with \(N_{S}\) labeled pixels (i.e., with its corresponding EMs and

Figure 1: Overview of the proposed method. A self-supervised learning module is used to generate supervised and unsupervised training data, \(\mathcal{D}_{U}\) and \(D_{S}\), to train the model. A generative model (with parameters \(\theta\)) is constructed with a prior model \(p_{\theta}(\mathbf{M},\mathbf{a},\mathbf{Z})\), describing the a priori behavior of the EMs, abundances, and variability parameters, and a mixing model \(p_{\theta}(\mathbf{y}|\mathbf{M},\mathbf{a},\mathbf{Z})\), containing both a linear and a nonlinear component. To solve HU, an approximate posterior distribution (with parameters \(\phi\)) \(q_{\phi}(\mathbf{M},\mathbf{a},\mathbf{Z}|\mathbf{y})\) is constructed, and independence assumptions between the abundances and the EM variability parameters are used to disentangle these variables. A semi-supervised learning-based optimization objective is used to learn the parameters of both the generative model and of the approximate posterior by maximizing a lower bound on the log-likelihood of both the supervised and unsupervised data, in \(\mathcal{L}_{\mathrm{Data}}(\theta,\phi,\mathcal{D})\), and also regularizations promoting abundance sparsity (\(\mathcal{L}_{\mathrm{Sparse}}(\theta,\phi,\mathcal{D})\)) and controlling the amount of nonlinearity in both the generative and inference models. The learned posterior then provides the HU solution as the EMs and the abundances.
abundances), which will be generated directly from \(\mathcal{D}_{U}\) based on a self-supervised learning strategy that is described in detail in Section IV-B. We also denote by \(\mathcal{D}=\mathcal{D}_{U}\bigcup\mathcal{D}_{S}\) the full dataset. Functions without accent (e.g., \(\mathbf{x}\)) belong to the mixing model, while functions with the _breve_ accent (e.g., \(\tilde{\mathbf{x}}\)) belong to the inference model. To simplify the notation, the pixel index will be omitted when possible.

_Section overview:_ In the following, we first describe the proposed mixing model in Section III-A. Then, the unmixing problem is formulated in a disentangled variational inference framework in Section III-B. The cost function for the semi-supervised learning objective is presented in III-C. Finally, the optimization of the proposed cost function is addressed in Section III-D. An overview of the proposed framework is given in Figure 1, which illustrates the interplay between the different parts of the statistical model representing the mixing (i.e., generative) and unmixing (i.e., inference) process, the generation of the training data, and the different terms in the loss function which is optimized to perform unmixing.

### _The mixing model_

_Abundance prior:_ When assigning a statistical distribution to the abundances, it is very important to account for the physical constraints. The Dirichlet distribution is a natural choice for modeling abundance vectors since it enforces the non-negativity and sum-to-one constraints of the model, thus, being appropriate to represent fractions. In this work we consider a flat Dirichlet distribution:

\[p(\mathbf{a})=\mathrm{Dir}(\mathbf{a};\mathbf{1}_{P})\,, \tag{2}\]

where \(\mathbf{1}_{P}\in\mathbb{R}^{P}\) is a \(P\)-dimensional vector of ones that contains its concentration parameters of the distribution. The flatness indicates the lack of prior knowledge over the abundances other than its physical constraints. The Dirichlet distribution properties made it a popular choice of prior for Bayesian HU strategies [82, 83, 84], and has also been used in conjunction with a Markov transition model to represent the abundances multitemporal HU [85].

_Endmember model:_ Several models have been proposed to account for EM variability by representing the spectral signatures of the materials in each pixel using deterministic or statistical models. Deterministic models include the use of additive [50] or multiplicative [51, 52, 86] perturbations of reference EM signatures, or a combination thereof [53]. The Gaussian [87], Beta [47] or Gaussian mixtures [48] distributions have been considered as statistical models. However, in practice EM variability can be very complex and specifying a distribution \(p(\mathbf{M})\) directly is very difficult. Moreover, these models do not explore an important property of spectral variability. Since the spectrum of most materials is a function of a small number of physico-chemical parameters [88, 89], the EM signatures are supported on a low-dimensional submanifold of the high-dimensional spectral space.

In this work, we consider a deep generative EM model to provide a flexible representation of EM signatures while accounting for their low intrinsic dimension [11]. Specifically, we model \(\mathbf{M}\) as a random variable that may be arbitrary but which follows a Gaussian distribution when conditioned on a set of low-dimensional latent variables \(\mathbf{Z}=[\mathbf{z}_{1},\dots,\mathbf{z}_{P}]\), \(\mathbf{z}_{k}\in\mathbb{R}^{H}\) (i.e., \(p(\mathbf{M}|\mathbf{Z})\) is Gaussian). The latent variables \(\mathbf{Z}\) control the variability of the spectral signatures, their dimension \(H\) being related to the flexibility of the model to represent the EMs in an HI. Results in [11] indicate that small values of \(H\) are sufficient for obtaining good unmixing results and represent EM variability well. In this way, although the conditional distribution is tractable, the marginal PDF \(p(\mathbf{M})=\int p(\mathbf{M}|\mathbf{Z})p(\mathbf{Z})d\mathbf{Z}\) can be arbitrary, rendering the model very flexible. Nonetheless, such a marginalization will not be necessary; instead, the latent variables \(\mathbf{Z}\) become an additional parameter that will be computed during inference. Considering the different EMs to be conditionally independent, the model becomes:

\[p_{\theta}(\mathbf{M}|\mathbf{Z})=\prod_{k=1}^{P}p_{\theta}(\mathbf{m}_{k}|\mathbf{z}_{k})\,, \tag{3}\]

with \(\mathbf{m}_{k}\in\mathbb{R}^{L}\) being the \(k\)-th column of \(\mathbf{M}\), and

\[p_{\theta}(\mathbf{m}_{k}|\mathbf{z}_{k})=\mathcal{N}\big{(}\mathbf{m}_{k};\mathbf{\mu}_{ \theta}^{m,k}(\mathbf{z}_{k}),\mathrm{diag}(\mathbf{\sigma}_{\theta}^{m,k}(\mathbf{z}_{k}) )\big{)}\,. \tag{4}\]

Here \(\mathcal{N}(\mathbf{x};\mathbf{\mu},\mathbf{\Sigma})\) denotes a Gaussian distribution of variable \(\mathbf{x}\) with mean \(\mathbf{\mu}\) and covariance matrix \(\mathbf{\Sigma}\). Functions \(\mathbf{\mu}_{\theta}^{m,k}(\mathbf{z}_{k})\) and \(\mathrm{diag}(\mathbf{\sigma}_{\theta}^{m,k}(\mathbf{z}_{k}))\) return the mean and (diagonal) covariance matrix of the distribution, and \(\theta\) is the set of parameters of the generative model. We assign a prior for \(\mathbf{Z}\) as

\[p(\mathbf{Z})=\prod_{k=1}^{P}p(\mathbf{z}_{k})\,,\quad p(\mathbf{z}_{k})=\mathcal{N}(\mathbf{z} _{k};\mathbf{0},\mathbf{I})\,. \tag{5}\]

Note that the number of EMs, \(P\), is assumed to be known a priori. In practice, it can be estimated from the HI using virtual dimensionality estimation techniques [89]. However, special care has to be taken to avoid overestimating the number of EMs in the scene, as explained in [90].

_Mixing model:_ Several models have been proposed to represent both macroscopic and intimate nonlinear mixtures, such as the BLMM and Hapke models [38, 4]. However, due to the complexity of the mixing process, specifying a precise model can be difficult. This has motivated the development of HU method based on nonparametric models where the nonlinearity is learned from the data, such as in kernel-based methods [9, 41] and nonlinear AEC networks [10, 16, 17]. In this work, we consider a decomposition of the nonlinear mixing process as the sum of a linear term (the LMM) and a nonparametric nonlinear component [9, 17, 91, 92]:

\[\mathbf{y}=\mathbf{M}\mathbf{a}+\mathbf{\mu}_{\theta}^{y}(\mathbf{a},\mathbf{M})+\mathbf{e}\,, \tag{6}\]

where \(\mathbf{\mu}_{\theta}^{y}\) denotes the nonlinear contribution. This model aims at representing both macroscopic and intimate mixtures [91], with the linear term representing the main contribution of macroscopic interactions, and the nonlinear term representing the effects of multiple scattering and intimate mixtures. Its nonparametric form does not make direct assumptions about the type of mixture in an HI, at the cost of having to be learned. Moreover, the amount of nonlinearity can be controlled by penalizing the norm of \(\mathbf{\mu}_{\theta}^{y}\) during the learning process, which becomes close to the LMM when \(\mathbf{\mu}_{\theta}^{y}\) is small.
Considering the noise \(\mathbf{e}\in\mathbb{R}^{L}\) to be independent, white and Gaussian, the data likelihood becomes:

\[p_{\theta}(\mathbf{y}|\mathbf{a},\mathbf{M})=\mathcal{N}\big{(}\mathbf{y};\mathbf{M}\mathbf{a}+\mathbf{\mu} _{\theta}^{y}(\mathbf{a},\mathbf{M}),\sigma_{\theta}^{y}\mathbf{I}\big{)}\,, \tag{7}\]

where \(\sigma_{\theta}^{y}\in\mathbb{R}_{+}\) is the noise variance in each band.

Since we assume that the abundances, the EMs and the noise are independent, this leads to the following factorization for the joint distribution:

\[p_{\theta}(\mathbf{y},\mathbf{a},\mathbf{M},\mathbf{Z})=p_{\theta}(\mathbf{y}|\mathbf{a},\mathbf{M})p_{ \theta}(\mathbf{a})p_{\theta}(\mathbf{M}|\mathbf{Z})p(\mathbf{Z})\,. \tag{8}\]

### _The unmixing problem_

The unmixing problem consists of finding the posterior distribution \(p(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})\). However, due to the choice of distributions in (8) computing the posterior analytically is intractable, and some approximations have to be performed. To this end, we use a variational approximation, which consists of specifying a parametric distribution \(q_{\phi}(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})\) from a sufficiently rich family, and finding the set of parameters \(\phi\) which makes it as close as possible to the true posterior distribution \(p(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})\)[93]. We consider the following factorization for this approximation:

\[q_{\phi}(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})=q_{\phi}(\mathbf{a}|\mathbf{y},\mathbf{M},\mathbf{Z})q_ {\phi}(\mathbf{M}|\mathbf{Z},\mathbf{y})q_{\phi}(\mathbf{Z}|\mathbf{y}). \tag{9}\]

Moreover, we also consider disentanglement learning within a statistical framework, simplifying the model by separating latent variables through statistical independence assumptions [20]. Thus, we assume that \(\mathbf{a}\) is independent of \(\mathbf{Z}\) conditioned on \(\mathbf{y}\) and \(\mathbf{M}\), and that \(\mathbf{M}\) is independent of \(\mathbf{y}\) conditioned on \(\mathbf{Z}\), that is:

\[q_{\phi}(\mathbf{a}|\mathbf{y},\mathbf{M},\mathbf{Z}) =q_{\phi}(\mathbf{a}|\mathbf{y},\mathbf{M})\,, \tag{10}\] \[q_{\phi}(\mathbf{M}|\mathbf{Z},\mathbf{y}) =q_{\phi}(\mathbf{M}|\mathbf{Z})\,. \tag{11}\]

Through these two assumptions, \(\mathbf{a}\) is disentangled from \(\mathbf{Z}\), and \(\mathbf{M}\) is disentangled from \(\mathbf{y}\). In the following, we define specific forms for each of those conditional PDFs.

For \(q_{\phi}(\mathbf{a}|\mathbf{y},\mathbf{M})\), we consider another Dirichlet distribution:

\[q_{\phi}(\mathbf{a}|\mathbf{y},\mathbf{M})=\mathrm{Dir}(\mathbf{a};\tilde{\mathbf{\gamma}}_{\phi} ^{a}(\mathbf{y},\mathbf{M}))\,, \tag{12}\]

where function \(\tilde{\mathbf{\gamma}}_{\phi}^{a}\) computes its concentration parameters.

Most works (such as AECs) consider "black-box" models to compute the abundances from the mixed pixels. However, it is important to incorporate prior knowledge about the model to improve the interpretability of the results when parametrizing the abundances posterior distribution in (12). Specifically, it has been shown that for autoencoder networks, the amount of nonlinearity in the mixing model and in the inference model are closely connected [17] in such a way that it is helpful to also decompose the encoder into a linear and nonlinear parts, so that the amount of nonlinearity in both parts of the model can be controlled. Thus, we propose to write \(\tilde{\mathbf{\gamma}}_{\phi}^{a}(\mathbf{y},\mathbf{M})\) as:

\[\tilde{\mathbf{\gamma}}_{\phi}^{a}(\mathbf{y},\mathbf{M})=s_{\text{ReLU}}\left(\tilde{ \mathbf{\gamma}}_{\phi}^{a,\text{lin}}(\mathbf{y},\mathbf{M})+\tilde{\mathbf{\gamma}}_{\phi}^ {a,\text{lin}}(\mathbf{y},\mathbf{M})\right)\,, \tag{13}\]

where \(\tilde{\mathbf{\gamma}}_{\phi}^{a,\text{lin}}(\mathbf{y},\mathbf{M})\) is a (piecewise) linear function which estimates the abundance mean and uncertainty based on a parametric model architecture which can also promote abundance sparsity effectively (as will be discussed in the following), while \(\tilde{\mathbf{\gamma}}_{\phi}^{a,\text{lin}}(\mathbf{y},\mathbf{M})\) is a nonparametric function that is able to compensate other nonlinearities in the model. Function \(s_{\text{ReLU}}(\mathbf{x})=\max(\mathbf{0},\mathbf{x})\) is the ReLU activation, which is used to ensure the nonnegativity of \(\tilde{\mathbf{\gamma}}_{\phi}^{a}(\mathbf{y},\mathbf{M})\). Thus, we can use regularizations to explicitly control the contribution of the nonparametric term \(\tilde{\mathbf{\gamma}}_{\phi}^{a,\text{lin}}(\mathbf{y},\mathbf{M})\) in the abundance estimates.

Distribution \(q_{\phi}(\mathbf{Z}|\mathbf{y})\), is assumed to factorize as \(q_{\phi}(\mathbf{Z}|\mathbf{y})=\prod_{k=1}^{P}q_{\phi}(\mathbf{z}_{k}|\mathbf{y})\), where each \(q_{\phi}(\mathbf{z}_{k}|\mathbf{y})\) follows a Gaussian distribution:

\[q_{\phi}(\mathbf{z}_{k}|\mathbf{y})=\mathcal{N}\big{(}\mathbf{z}_{k};\tilde{\mathbf{\mu}}_{ \phi}^{z,k}(\mathbf{y}),\mathrm{diag}(\tilde{\mathbf{\sigma}}_{\phi}^{z,k}(\mathbf{y})) \big{)}\,, \tag{14}\]

in which functions \(\tilde{\mathbf{\mu}}_{\phi}^{z,k}\) and \(\tilde{\mathbf{\sigma}}_{\phi}^{z,k}\) compute its mean and the elements of its diagonal covariance matrix, respectively.

We also assume that \(q_{\phi}(\mathbf{M}|\mathbf{Z})\) can be factorized as \(q_{\phi}(\mathbf{M}|\mathbf{Z})=\prod_{k=1}^{P}q_{\phi}(\mathbf{m}_{k}|\mathbf{z}_{k})\), with:

\[q_{\phi}(\mathbf{m}_{k}|\mathbf{z}_{k})=\mathcal{N}\big{(}\mathbf{m}_{k};\tilde{\mathbf{\mu}}_{ \phi}^{m,k}(\mathbf{z}_{k}),\mathrm{diag}(\tilde{\mathbf{\sigma}}_{\phi}^{m,k}(\mathbf{z}_ {k}))\big{)}\,. \tag{15}\]

#### Sharing parameters

To simplify the inference process, we consider the same form for the EM conditional distribution in both (4) and (15). More precisely, we use \(\tilde{\mathbf{\mu}}_{\phi}^{m,k}(\mathbf{z}_{k})=\mathbf{\mu}_{\theta}^{m,k}(\mathbf{z}_{k})\) and \(\tilde{\mathbf{\sigma}}_{\phi}^{m,k}(\mathbf{z}_{k})=\mathbf{\sigma}_{\theta}^{m,k}(\mathbf{z}_ {k})\), therefore making \(q_{\phi}(\mathbf{M}|\mathbf{Z})=p_{\theta}(\mathbf{M}|\mathbf{Z})\). This is also known as "bottom-up" inference [94], and allows for parameter sharing to reduce the complexity of the model.

An illustration of the relationship between the different random variables, \(\mathbf{a},\mathbf{Z},\mathbf{M}\) and \(\mathbf{y}\), in both the generative and inference models according to the selected factorizations of the PDFs and their parametrizations can be seen in Figure 2.

### _Objective function_

To learn the parameters of the model and perform HU, we maximize a cost function composed of three terms, defined as:

\[\mathcal{L}_{\mathrm{T}}(\theta,\phi;\mathcal{D})=\mathcal{L}_{\mathrm{Data}}( \theta,\phi;\mathcal{D})-\mathcal{L}_{\mathrm{Sparse}}(\theta,\phi;\mathcal{D})- \mathcal{R}(\theta,\phi)\,. \tag{16}\]

The first term attempts to fit the model to the semi-supervised training dataset \(\mathcal{D}\), and aims to maximize the data likelihood \(\log p(\mathcal{D})\). The second term aims to promote sparse abundance estimates, while the last term is a data-independent regularization on the model parameters.

Figure 2: Overview of the dependencies between the variables in the proposed generative and inference (unmixing) models. Boxes denote functions, while the blue and red lines/arrows denote the data generating model and the inference model, respectively. _Data generating model:_ A latent variable \(\mathbf{Z}\) generates the EM matrix \(\mathbf{M}\), which together with the abundances \(\mathbf{a}\) generates the mixed pixel by an additive linear/nonlinear model. _Inference model:_ The learned model maps from the observed pixels to the latent EM representation \(\mathbf{Z}\), and from these to the estimated abundances using a “bottom-up” inference strategy and a two-stream neural network, respectively.
#### Iii-C1 **Term \(\mathcal{L}_{\rm Data}(\theta,\phi;\mathcal{D})\)**

The main criterion used to learn the parameters of the model is to maximize the likelihood of the data:

\[\log p(\mathcal{D})=\sum_{\mathbf{y}\in\mathcal{D}_{U}}\log p(\mathbf{y})+ \sum_{(\mathbf{y},\mathbf{a},\mathbf{M})\in\mathcal{D}_{S}}\log p(\mathbf{y},\mathbf{a},\mathbf{M})\,, \tag{17}\]

where the first term in the right hand side of (17) consists in the data likelihood for the unlabeled dataset \(\mathcal{D}_{U}\), while the second term is the data likelihood for the labeled dataset \(\mathcal{D}_{S}\) (i.e., in which variables \(\mathbf{a}\) and \(\mathbf{M}\) are treated as fixed/observed, and thus their likelihood is also maximized).

However, maximizing the expression in (17) is intractable because it involves marginalizing different variables in the joint distribution (8), which is a common issue in latent variable models [95]. This problem can be circumvented by maximizing the so-called _evidence lower bound_ (ELBO) of each term in (17):

\[\sum_{\mathbf{y}\in\mathcal{D}_{U}}\log p(\mathbf{y})\geq L_{U}(\theta, \phi;\mathcal{D}_{U})\,, \tag{18}\] \[\sum_{(\mathbf{y},\mathbf{a},\mathbf{M})\in\mathcal{D}_{S}}\log p(\mathbf{y},\bm {a},\mathbf{M})\geq L_{S}(\theta,\phi;\mathcal{D}_{S})\,, \tag{19}\]

where \(L_{U}(\theta,\phi;\mathcal{D}_{U})\) and \(L_{S}(\theta,\phi;\mathcal{D}_{S})\) are given by [20, 25]:

\[L_{U}(\theta,\phi; \mathcal{D}_{U})=\sum_{\mathbf{y}\in\mathcal{D}_{U}}\mathbb{E}_{q_{ \phi}(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})}\Big{[}\log p_{\theta}(\mathbf{y}|\mathbf{a},\mathbf{M}, \mathbf{Z})\] \[+\log p_{\theta}(\mathbf{Z},\mathbf{a},\mathbf{M})-\log q_{\phi}(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})\Big{]}\,, \tag{20}\]

and

\[L_{S}(\theta,\phi; \mathcal{D}_{S})=\sum_{(\mathbf{y},\mathbf{a},\mathbf{M})\in\mathcal{D}_{S}} \mathbb{E}_{q_{\phi}(\mathbf{Z}|\mathbf{y},\mathbf{a},\mathbf{M})}\Big{[}\log p_{\theta}(\mathbf{ y}|\mathbf{a},\mathbf{M},\mathbf{Z})\] \[+\log p_{\theta}(\mathbf{Z},\mathbf{a},\mathbf{M})-\log q_{\phi}(\mathbf{Z}|\mathbf{ y},\mathbf{a},\mathbf{M})\Big{]}\,. \tag{21}\]

Moreover, it can be shown that maximizing the ELBO equivalently minimizes the KL divergence between the variational posterior \(q_{\phi}(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})\) and the (unknown) true posterior \(p_{\theta}(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})\)[96].

The cost function to be maximized is a combination of the lower bounds \(L_{U}\) and \(L_{S}\), and an additional regularization term [20, 25]:

\[\mathcal{L}_{\rm Data}(\theta,\phi;\mathcal{D}) =L_{U}(\theta,\phi;\mathcal{D}_{U})+\lambda\Big{[}L_{S}(\theta, \phi;\mathcal{D}_{S})\] \[+\beta\sum_{(\mathbf{y},\mathbf{a},\mathbf{M})\in\mathcal{D}_{S}}\log q_{ \phi}(\mathbf{a},\mathbf{M}|\mathbf{y})\Big{]}\,, \tag{22}\]

where parameter \(\lambda\in\mathbb{R}_{+}\) is used to balance the contribution of the supervised and unsupervised terms in the cost function. This is especially important when the size of the supervised dataset is much smaller than that of the unsupervised one [20]. The additional term \(\log q_{\phi}(\mathbf{a},\mathbf{M}|\mathbf{y})\), weighted by \(\beta\in\mathbb{R}_{+}\), is a regularization added to the cost function to promote models in which the labeled data (abundances and EMs) have high likelihood in the posterior, conditioned on the pixels [25].

#### Iii-C2 **Term \(\mathcal{L}_{\rm Sparse}(\theta,\phi;\mathcal{D})\)**

An important property of many HIs is that while many materials may be present in the whole scene, only a small subset of them are usually present in each pixel. This leads the resulting abundance vectors to be sparse, with only a few of its elements being significantly far from zero [29]. This property can be used to improve the conditioning of the unmixing problem in the form of regularization strategies. In this work, we penalize the concentration parameters of the abundance posterior distribution (12) to encourage sparse abundance reconstructions:

\[\mathcal{L}_{\rm Sparse}(\theta,\phi;\mathcal{D})=\tau\sum_{(\bm {y},\mathbf{a},\mathbf{M})\in\mathcal{D}_{S}}\left\|\tilde{\gamma}_{\phi}^{a}(\mathbf{y}, \mathbf{M})\right\|_{1/2}\] \[+\tau\sum_{\mathbf{y}\in\mathcal{D}_{U}}\mathbb{E}_{q_{\phi}(\mathbf{M}| \mathbf{Z})q_{\phi}(\mathbf{Z}|\mathbf{y})}\Big{\{}\big{\|}\tilde{\gamma}_{\phi}^{a}(\mathbf{y},\mathbf{M})\|_{1/2}\Big{\}}\,, \tag{23}\]

where the averages of the L\({}_{1/2}\) semi-norm of the Dirichlet concentration parameters are performed with respect to the supervised (first term) and unsupervised (second term) data portions, and \(\tau\in\mathbb{R}_{+}\) is a parameter controlling the regularization. Furthermore, we considered the L\({}_{1/2}\) norm penalization as it is a more effective measure of sparseness when compared to the traditionally used L\({}_{1}\) norm penalty [24].

#### Iii-C3 **Term \(\mathcal{R}(\theta,\phi)\)**

To control the nonlinear contributions and retain the interpretability of both the encoder (6) and decoder (13) models, we add as a penalization term

\[\mathcal{R}(\theta,\phi)=\varsigma_{1}\big{\|}\mathbf{\mu}_{\theta}^{y}\big{\|}_{ \rm FNN}+\varsigma_{2}\big{|}\tilde{\gamma}_{\phi}^{a,\rm lin}\big{\|}_{\rm FNN}\,, \tag{24}\]

where \(\|f\|_{\rm FNN}\) denotes the norm of some feedforward neural network \(f\), defined as the sum of the norms of its parameters at all layers. For instance, if \(f\) represents a multilayer perceptron (MLP) with weights \(\mathbf{W}_{\ell}\) and biases \(\mathbf{b}_{\ell}\) at layer \(\ell\), the norm is computed as \(\|f\|_{\rm FNN}=\sum_{\ell}|\mathbf{W}_{\ell}\|_{F}+\|\mathbf{b}_{\ell}\|_{2}\), with \(\|\cdot\|_{F}\) denoting the Frobenius norm. The first regularization term in (24) penalizes the amount of nonlinearity in the mixing model (6), while the second term penalizes the contribution of the nonparametric component in the abundance estimates in (13). Constants \(\varsigma_{1},\varsigma_{2}\in\mathbb{R}_{+}\) are regularization parameters.

### _Optimizing the cost function_

Due to the nature of the probabilistic model considered in this work, directly maximizing the cost function \(\mathcal{L}_{\rm T}(\theta,\phi;\mathcal{D})\) in (16) is not trivial. This occurs because it involves some conditional distributions which are not available in the model factorization devised in Section III-B. Moreover, it also requires optimizing expectations with respect to both Gaussian and Dirichlet distributions depending on the model parameters. We address these problems in the following.

#### Iii-D1 **Rewriting the unavailable distributions**

The first difficulty comes from the cost function \(\mathcal{L}_{\rm Data}(\theta,\phi;\mathcal{D})\) in (22). First, it requires sampling from and evaluating \(q_{\phi}(\mathbf{Z}|\mathbf{y},\mathbf{a},\mathbf{M})\) (see (21)), which is not available under the selected posterior factorization in (9). Similarly, we also do not have access to the marginal posterior \(q_{\phi}(\mathbf{a},\mathbf{M}|\mathbf{y})\). To circumvent the first issue,
we can proceed as in [20] and rewrite the last two terms of \(\mathcal{L}_{\mathrm{Data}}(\theta,\phi;\mathcal{D})\) as:

\[L_{S}(\theta,\phi;\mathcal{D}_{S})+\beta\sum_{(\mathbf{y},\mathbf{a}, \mathbf{M})\in\mathcal{D}_{S}}\log q_{\phi}(\mathbf{a},\mathbf{M}|\mathbf{y})\] \[=\sum_{(\mathbf{y},\mathbf{a},\mathbf{M})\in\mathcal{D}_{S}}\Big{[}\mathbb{E}_ {\mathbb{E}_{q_{\phi}}(\mathbf{Z}|\mathbf{y},\mathbf{a},\mathbf{M})}\Big{\{}\log p_{\theta}( \mathbf{y}|\mathbf{a},\mathbf{M},\mathbf{Z})\] \[+\log p_{\theta}(\mathbf{a},\mathbf{M},\mathbf{Z})-\log q_{\phi}(\mathbf{a},\mathbf{ M},\mathbf{Z}|\mathbf{y})\Big{\}}\] \[+(1+\beta)\log q_{\phi}(\mathbf{a},\mathbf{M}|\mathbf{y})\Big{]}\,, \tag{25}\]

where we added and subtracted the term \(\log q_{\phi}(\mathbf{a},\mathbf{M}|\mathbf{y})\). Note that the variational posterior inside the expectation is now in the same form as (9).

However, this still requires optimizing an expectation taken with respect to the unavailable \(q_{\phi}(\mathbf{Z}|\mathbf{y},\mathbf{a},\mathbf{M})\). This can be approximated using importance sampling [20], which allows us to approximate the expectation w.r.t. this term by sampling from \(\mathbf{Z}^{(i)}\sim q_{\phi}(\mathbf{Z}|\mathbf{y})\) (for each datapoint \(\mathbf{y}\in\mathcal{D}_{U}\)), which is available under the factorization (9). Using self-normalized importance sampling with \(K\) samples, we can write the expectation in the right hand side of (25) as:

\[\mathbb{E}_{q_{\phi}(\mathbf{Z}|\mathbf{y},\mathbf{a},\mathbf{M})}\big{\{}\log p _{\theta}(\mathbf{y},\mathbf{a},\mathbf{M},\mathbf{Z})-\log q_{\phi}(\mathbf{Z},\mathbf{a},\mathbf{M}|\mathbf{ y})\big{\}}\] \[=\mathbb{E}_{q_{\phi}(\mathbf{Z}|\mathbf{y})}\bigg{\{}\frac{q_{\phi}(\bm {Z}|\mathbf{y},\mathbf{a},\mathbf{M})}{q_{\phi}(\mathbf{Z}|\mathbf{y})}\bigg{[}\log\frac{p_{\theta} (\mathbf{y},\mathbf{a},\mathbf{M},\mathbf{Z})}{q_{\phi}(\mathbf{Z},\mathbf{a},\mathbf{M}|\mathbf{y})}\bigg{]} \bigg{\}}\] \[=\mathbb{E}_{q_{\phi}(\mathbf{Z}|\mathbf{y})}\bigg{\{}\frac{q_{\phi}(\bm {Z}|\mathbf{y},\mathbf{a},\mathbf{M})q_{\phi}(\mathbf{M},\mathbf{a}|\mathbf{y})}{q_{\phi}(\mathbf{Z}|\mathbf{y} )q_{\phi}(\mathbf{M},\mathbf{a}|\mathbf{y})}\] \[\times\bigg{[}\log\frac{p_{\theta}(\mathbf{y},\mathbf{a},\mathbf{M},\mathbf{Z})}{ q_{\phi}(\mathbf{Z},\mathbf{a},\mathbf{M}|\mathbf{y})}\bigg{]}\bigg{\}}\] \[=\mathbb{E}_{q_{\phi}(\mathbf{Z}|\mathbf{y})}\bigg{\{}\frac{q_{\phi}(\bm {M}|\mathbf{Z})}{q_{\phi}(\mathbf{M}|\mathbf{y})}\bigg{[}\log\frac{p_{\theta}(\mathbf{y},\mathbf{ a},\mathbf{M},\mathbf{Z})}{q_{\phi}(\mathbf{Z},\mathbf{a},\mathbf{M}|\mathbf{y})}\bigg{]}\bigg{\}}\] \[\simeq\sum_{i=1}^{K}\frac{\omega^{(i)}}{\Omega}\log\frac{p_{\theta }(\mathbf{y},\mathbf{a},\mathbf{M},\mathbf{Z}^{(i)}}{q_{\phi}(\mathbf{Z}^{(i)},\mathbf{a},\mathbf{M}|\mathbf{y })}\,, \tag{26}\]

where (10) and (11) were used to obtain the third equality. Variable \(\mathbf{Z}^{(i)}\) is sampled as \(\mathbf{Z}^{(i)}\sim q_{\phi}(\mathbf{Z}|\mathbf{y})\), and the importance weights \(\omega^{(i)}\) and their normalization factor \(\Omega\) are obtained as

\[\frac{q_{\phi}(\mathbf{M}|\mathbf{Z}^{(i)})}{q_{\phi}(\mathbf{M}|\mathbf{y})} \varpropto q_{\phi}(\mathbf{M}|\mathbf{Z}^{(i)})\triangleq\omega^{(i)}\,, \tag{27}\] \[\Omega =\sum_{i=1}^{K}\omega^{(i)}\,. \tag{28}\]

Note that since \(q_{\phi}(\mathbf{M}|\mathbf{y})\) is constant in (27), i.e., it does not depend on \(i\), and is also present in \(\Omega\) in (26) it cancels in (28). Therefore, \(q_{\phi}(\mathbf{M}|\mathbf{y})\) (which is not available under the chosen posterior factorization) is not required.

Finally, to approximate term \(\log q_{\phi}(\mathbf{a},\mathbf{M}|\mathbf{y})\) in (25), note that it can be lower bounded as [20]:

\[\log q_{\phi}(\mathbf{a}, \mathbf{M}|\mathbf{y})\geq\mathbb{E}_{q_{\phi}(\mathbf{Z}|\mathbf{y})}\bigg{\{} \log\frac{q_{\phi}(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})}{q_{\phi}(\mathbf{Z}|\mathbf{y})} \bigg{\}} \tag{29}\] \[=\log q_{\phi}(\mathbf{a}|\mathbf{M},\mathbf{y})+\mathbb{E}_{q_{\phi}(\mathbf{Z}| \mathbf{y})}\Big{\{}\log q_{\phi}(\mathbf{M}|\mathbf{Z})\Big{\}}\,.\]

Then, using a Monte Carlo estimator of the expectation in (29) with \(K\) samples gives us:

\[\mathbb{E}_{q_{\phi}(\mathbf{Z}|\mathbf{y})}\Big{\{}\log q_{\phi}(\mathbf{M}| \mathbf{Z})\Big{\}} \simeq\frac{1}{K}\sum_{i=1}^{K}\log q_{\phi}(\mathbf{M}|\mathbf{Z}^{(i)})\] \[=\frac{1}{K}\sum_{i=1}^{K}\log\omega^{(i)}\,. \tag{30}\]

It can be seen that this approximation can be computed directly using the importance weights obtained in (27).

To achieve the final cost function we re-write the terms in (16) leveraging the previous results for the different terms. For the terms in \(\mathcal{L}_{\mathrm{Data}}(\theta,\phi;\mathcal{D})\) shown in (25) we used the lower bound in (29) and the importance sampling Monte Carlo approximation in (26). For the remaining \(L_{U}(\theta,\phi;\mathcal{D}_{U})\) term and \(\mathcal{L}_{\mathrm{Sparse}}(\theta,\phi;\mathcal{D})\) we applied Monte Carlo approximations with \(K_{E}\) samples to compute its expectations. The resulting cost function is shown in equation (31).

### _Reparametrization_

In order to optimize the cost function in (31) with respect to the parameters \(\theta\) and \(\phi\), it is necessary to estimate the gradient of expressions involving variables \(\mathbf{Z}^{(i)}\), \(\mathbf{M}^{(i)}\) and \(\mathbf{a}^{(i)}\), which are sampled from the variational posterior distribution \(q_{\phi}(\mathbf{a},\mathbf{M},\mathbf{Z}|\mathbf{y})\). However, since this distribution depends on the parameters \(\phi\), it is necessary to account for the dependency of the samples \(\mathbf{Z}^{(i)}\), \(\mathbf{M}^{(i)}\) and \(\mathbf{a}^{(i)}\) on \(\phi\) when computing the gradients. To this end, we consider the reparametrization trick, which provides low-variance gradient estimates for such problems [97]. In general, this is performed by writing the considered random variables (say, \(\mathbf{x}\)) in terms of a distribution that does not depend on \(\phi\), e.g., for a given function \(f\) to be differentiated, we have:

\[f(\mathbf{x})=f(g(\mathbf{\epsilon},\phi))\,, \tag{32}\]

where \(g\) is a function such that \(x\) and \(g(\mathbf{\epsilon},\phi)\) have the same distribution, and \(\mathbf{\epsilon}\sim p(\mathbf{\epsilon})\) does not depend on \(\phi\). Thus, the gradient \(\partial f(\mathbf{x})/\partial\phi\) can be computed straightforwardly and estimated using Monte Carlo sampling.

This strategy can be applied to reparametrize the samples from \(q_{\phi}(\mathbf{M}|\mathbf{Z})\) and \(q_{\phi}(\mathbf{Z}|\mathbf{y})\). Since the variables \(\mathbf{Z}\) and \(\mathbf{M}\) are conditionally Gaussian, the reparametrization trick is straightforward to apply (as in [98]), resulting in:

\[\mathbf{Z}^{(i)} =\big{[}\tilde{\mathbf{\mu}}_{\phi}^{z,1}(\mathbf{y}),\ldots,\tilde{\mathbf{\mu }}_{\phi}^{z,P}(\mathbf{y})\big{]}\big{\circ}\,\mathbf{\Xi}_{z}^{(i)}\,, \tag{33}\] \[\mathbf{M}^{(i)} =\big{[}\tilde{\mathbf{\mu}}_{\phi}^{m,1}(\mathbf{Z}_{1}^{(i)}),\ldots, \tilde{\mathbf{\mu}}_{\phi}^{m,P}(\mathbf{Z}_{P}^{(i)})\big{]}\] \[\quad+\big{[}\tilde{\mathbf{\sigma}}_{\phi}^{m,1}(\mathbf{Z}_{1}^{(i)}), \ldots,\tilde{\mathbf{\sigma}}_{\phi}^{m,P}(\mathbf{Z}_{P}^{(i)})\big{]}\big{\circ}\, \mathbf{\Xi}_{m}^{(i)}\,, \tag{34}\]

where \(\circ\) is the Hadamard (i.e., elementwise) product and \(\mathbf{\
in the softmax basis [99] or reparametrizing rejection samplers [22]. Here we consider the so-called pathwise derivative estimator proposed in [23], which allows us to estimate the derivative of a stochastic variable \(\mathbf{a}\) sampled from a Dirichlet distribution \(\mathrm{Dir}(\mathbf{a};\mathbf{\gamma})\) as:

\[\frac{\partial a_{i}}{\partial\gamma_{j}}=-\frac{\frac{\partial E_{ \mathrm{Beta}}}{\partial\gamma_{j}}\big{(}a_{j};\gamma_{j},\sum_{k=1}^{P} \gamma_{k}-\gamma_{j}\big{)}}{f_{\mathrm{Beta}}\big{(}a_{j};\gamma_{j},\sum_{k =1}^{P}\gamma_{k}-\gamma_{j}\big{)}}\Big{(}\frac{\delta_{i-j}-a_{i}}{1-a_{j}} \Big{)}\,, \tag{35}\]

where \(f_{\mathrm{Beta}}\) and \(F_{\mathrm{Beta}}\) are the Beta probability and cumulative distribution functions, respectively, and \(\delta_{i-j}\) the Kronecker delta function, satisfying \(\delta_{i-j}=1\) if \(i=j\) and zero otherwise (see [23] for more details). This expression can be used to compute the gradient of the cost function using the chain rule.

### _On the interpretability of the proposed approach_

Much interest has been recently dedicated to interpretability in the field of machine learning [13, 100], and while a mathematical formalization of interpretability is difficult, a definition that is appropriate in this context is when it allows a user/human observer to understand the model and algorithm's behavior [100, 10]. In the proposed method, interpretability mainly comes from the statistical modeling and residual network design.

First, the unmixing problem tackled in this work involves multiple nonidealities, such as nonlinearity and variability of the endmembers. By modeling both the generative and inference (unmixing) models as PDFs (in (8) and (9)) and employing adequate conditional independence assumptions on the different model variables (e.g., in (10) (11)), we are able to separately account for each of these effects, instead of having a single black box model. Moreover, the independence assumptions used in the inference model allow us to disentangle the influence of each of these effects in the recovered abundances [20]. This framework also leads to a principled training procedure based on semi-supervised learning, where the contribution of the data with and without labels (which might be extracted from the HI through self-supervised learning, as detailed later) is clear. Another point is that the residual formulation of the mixing model, i.e., a linear mixing model augmented with a flexible nonlinear component, see (6), copes with different types of mixtures, such as, e.g., the multiple mixture pixel model [91] or BLMMs [4].

## IV Neural network design and training

In this section, we detail different aspects of the method, consisting in: 1) the choices of model and NN architectures, 2) the construction of the supervised training data using a self-supervision strategy, and 3) the cost function optimization strategy.

### _Model details and neural network architectures_

#### Iv-A1 **Endmember networks**

For the EM networks, which is used in both the distribution \(p_{\theta}(\mathbf{m}_{k}|\mathbf{z}_{k})\) in the mixing model (4) and \(q_{\phi}(\mathbf{m}_{k}|\mathbf{z}_{k})\) in the inference network (15), we considered the following architectures. For the mean vectors \(\tilde{\mathbf{\mu}}_{\phi}^{m,k}(\mathbf{z}_{k})=\mathbf{\mu}_{\theta}^{m,k}(\mathbf{z}_{k})\), we considered a fully connected MLP as in [11] with six layers containing \(H\), \(\max\big{\{}[L/10],\,H+1\big{\}}\), \(\max\big{\{}[L/4],\,H+2\big{\}}+3\), \([1.2\times L]+5\) and \(L\) neurons. The ReLU activation function was used on the hidden layers and a sigmoid on the last layer. For the covariances, for simplicity we used an EM-wise learnable constant \(\sigma^{m,k}\in\mathbb{R}_{+}\) resulting in an isotropic model, which led to \(\tilde{\mathbf{\sigma}}_{\phi}^{m,k}(\mathbf{z}_{k})=\mathbf{\sigma}_{\theta}^{m,k}(\mathbf{z }_{k})=\sigma^{m,k}_{m}\mathbf{1}_{L}\), where \(\mathbf{1}_{L}\) is an \(L\)-dimensional vector of ones.

#### Iv-A2 **Variability inference network**

For the variational distributions \(q_{\phi}(\mathbf{z}_{k}|\mathbf{y})\) in (14), we considered two fully connected MLPs with partially shared parameters to compute the mean and diagonal covariance elements. For the mean \(\tilde{\mathbf{\mu}}_{\phi}^{z,k}(\mathbf{y})\), we considered four fully connected layers with \(L\), \(5H\), \(2H\) and \(H\) neurons. For the diagonal covariance elements \(\tilde{\mathbf{\sigma}}_{\phi}^{z,k}(\mathbf{y})\), we considered six layers with \(L\), \(5H\), \(2H\), \(2H\), \(2H\) and \(H\) neurons. The first two hidden layers were shared between both networks, and the ReLU activation function was used in all hidden layers, with a linear activation in the last layer.

#### Iv-A3 **Nonlinear mixing (decoder) network**

The measurement model \(p_{\theta}(\mathbf{y}|\mathbf{a},\mathbf{M})\) in (7) was parametrized as follows. For the function \(\mathbf{\mu}_{\theta}^{y}(\mathbf{a},\mathbf{M})\) representing the nonlinear contribution in the mixing model (6), we consider fully connected MLP as used in [17]. It consisted of five fully connected layers, with \(P(L+1)\), \(PL\), \(L\), \(L\) and \(L\) neurons. The ReLU activation was used in all hidden layers, and a linear activation was used in the last layer. The observation noise variance \(\sigma_{\theta}^{y}\) in (7) is set as a positive learnable constant.
#### Iv-A4 **Abundance posterior (encoder) network**

As discussed previously, we consider a two-stream architecture in (13) to estimate the abundance posterior distribution (12). For the piecewise linear reconstruction function \(\hat{\gamma}_{\phi}^{\alpha,\text{lin}}(\mathbf{y},\mathbf{M})\), we propose to consider an unrolled NN architecture that is suitable to sparse regression. Unrolling algorithms to create NNs consists of converting the iterative equations of traditional optimization algorithms (e.g., gradient descent or iterative shrinkage/thresholding) into NN layers. The different parameters of the algorithms, which are traditionally fixed (or come from a forward measurement model) become the trainable parameters of the NN [101]. Such architectures have become very prominent in machine learning solution to various inverse problems due to their high interpretability [101].

In this work, we consider an architecture inspired by the LISTA method [102] with \(M\) network layers, which unfolds the iterative shrinkage/thresholding algorithm for sparse regression. This leads to the following architecture, with the \(m\)-th hidden network layer, \(m\in\{1,\ldots,M-1\}\), given by:

\[\mathbf{h}^{(m+\frac{1}{2})} =\mathbf{h}^{(m)}-\eta^{(m)}\mathbf{M}^{\top}\big{(}\mathbf{M}\mathbf{h}^{(m)}- \mathbf{y}\big{)}\,, \tag{36}\] \[\mathbf{h}^{(m+1)} =s_{\text{ReLU}}\big{(}\mathbf{h}^{(m+\frac{1}{2})}-\eta^{\text{sp}} \eta^{(m)}\mathbf{1}_{P}\big{)}\,, \tag{37}\]

where \(\mathbf{h}^{(m)}\) denotes the signal in the \(m\)-th later, \(\eta^{(m)}\in\mathbb{R}_{+}\) is a learnable step size and \(\eta^{\text{sp}}\in\mathbb{R}_{+}\) a parameter that controls sparsity of the reconstruction. The hidden layers can be interpreted as a gradient step of a least squares problem, followed by a combined shrinkage and projection to the nonnegative orthant. The output of the first layer is computed as \(\mathbf{h}^{(1)}=\mathbf{M}^{\dagger}\mathbf{y}\), with \(\mathbf{M}^{\dagger}\) denoting the pseudoinverse of \(\mathbf{M}\). The last layer, computed as \(\mathbf{h}^{(M)}=\eta^{\text{unc}}\mathbf{h}^{(M-1)}\), consists of a simple scaling of the result by a scalar parameter \(\eta^{\text{unc}}\in\mathbb{R}_{+}\), which controls the spread of the Dirichlet distribution. This approach is interesting since the learnable parameters are just the step sizes \(\eta^{\text{unc}}\) and sparsity and uncertainty parameters \(\eta^{\text{sp}}\), \(\eta^{(m)}\), which maintains a high degree of interpretability in the architecture.

For the nonlinear part of the abundance posterior \(\hat{\gamma}_{\phi}^{\alpha,\text{lin}}(\mathbf{y},\mathbf{M})\), we considered an MLP with six fully connected layers containing \(L\), \(2L\), \(\lfloor L/2\rfloor\), \(\lfloor L/4\rfloor\), \(4P\) and \(P\) neurons, where \(\lfloor\cdot\rfloor\) denotes rounding to the nearest integer. The ReLU activation function was used in the hidden layers, and a linear activation in the last layer.

### _Self-supervised training strategy_

An important aspect of the semi-supervised framework is the need to obtain some amount of supervised training data. As done in previous works [10, 11], we propose to use a self-supervised strategy based on methods which extract different samples from EMs directly from an observed HI (i.e., in the form of pure pixels) [5, 103]. These methods have been originally proposed as image-based spectral library construction to address spectral variability in library-based HU algorithms without the requirement of laboratory or _in situ_ data collection. Nonetheless, these are also suitable strategies to generate training data for learning-based HU methods.

In this paper, we consider the strategy proposed in [11], to generate a supervised synthetic dataset \(\mathcal{D}_{\text{S}}\) directly from the HI being analyzed by following a two step procedure. In the first step, we create a dictionary \(\mathcal{D}_{\text{pps}}\) of pure pixels. Pure pixels are obtained by extracting a predetermined number \(N_{\text{pps}}\) of pixels from the HI which are spectrally close (e.g., having the smallest spectral angles) to a set of reference EMs obtained using a traditional EM extraction algorithm (such as the VCA [104]). The second step consists in an iterative algorithm for generating synthetic data that incorporates the variability existing within \(\mathcal{D}_{\text{pps}}\). For this, at each iteration \(k\) we generate tuples \(\{(\mathbf{y}_{j},\mathbf{a}_{j},\mathbf{M}_{k})\}_{j=1}^{P}\) where \(\mathbf{M}_{k}\in\mathbb{R}^{L\times P}\) is an EM matrix sampled from the pure pixel dictionary. For each \(k\) we generate \(P\) abundance vectors \(\mathbf{a}_{j}\), \(j\in\{1,\ldots,P\}\), with the elements \(a_{j,i}\) of the \(j\)-th vector satisfying \(a_{j,i}=1\) if \(i=p\) and \(a_{j,i}=0\) otherwise (i.e., a one-hot encoding of each EM). Then, the mixed pixels are generated according to \(\mathbf{y}_{j}=\mathbf{M}_{k}\mathbf{a}_{j}+\mathbf{e}\), with \(\mathbf{e}\) being white Gaussian noise.

Note that since the synthetic abundances are one-hot encodings, this leads to a simple but effective means of generating training data, since nonlinear interactions (i.e., a nonlinear mixing model) don't need to be specified a priori due to the pixels being pure pixels. Nevertheless, more general methods using a finer sampling of the unit simplex and different choices of mixing models can also be used to generate the training data at the expense of a higher amount of user supervision and possibly at a higher computation cost.

### _Optimization_

To optimize the cost function, we considered the Adam stochastic optimization method [26]. We used a batch size of \(16\) and trained the model for up to \(30\) epochs, stopping when the relative increase in the value of the cost function \(\hat{\mathcal{L}}_{\text{T}}(\theta,\phi;\mathcal{D})\) between two epochs was smaller than \(0.01\). The learning rate was initially set as \(0.001\) and multiplied by \(0.9\) after each epoch until the \(10\)th epoch using a learning rate scheduler, after which is was kept fixed.

## V Experiments

### _Experimental setup_

The proposed method is compared with the fully constrained least squares (FCLS) algorithm, with EMs extracted by the VCA [104], with the PLMM [50], ELMM [51], and GLMM [52] algorithms, which account for EM variability, with the NDU [105], which considers a nonparametric nonlinear mixing model, and with recent deep learning-based strategies that address EM variability and nonlinearity, being DeepGUn [11], RBF-AEC [69] and EGU-Net [10]. In each experiment the parameters of each algorithm were adjusted in order to obtain the best results. All experiments were run in an AMD Ryzen 5 3500U portable computer with 6GB of RAM.

The proposed method was implemented in Pytorch. We consider \(K_{E}=1\) sample for the Monte Carlo estimation of the expectations, and \(K=5\) for the importance sampling. The supervised dataset \(\mathcal{D}_{S}\) was constructed using the self-supervised strategy detailed in Section IV-B, with VCA being used to extract the reference EM signatures, \(N_{\text{pps}}=100\) pure pixels from each EM being extracted to generate the dataset, and a signal-to-noise ratio (SNR) of 30dB was considered for
the additive noise. The dimension of the EM latent space was set as \(H=2\), and \(M=11\) layers (i.e., ten hidden layers) were considered in the piecewise-linear abundance inference network detailed in Section IV-A4. The remaining parameters of the method were manually selected within the following intervals: \(\lambda\in[0.01,1000]\) with one value per decade; \(\tau\) was either zero or selected in \(\tau\in[0.001,1]\) with two values per decade; and the nonlinearity regularizations in the interval \(\varsigma_{1},\varsigma_{2}\in[10^{-4},10^{5}]\) with one value per decade. For all methods, the number of EMs \(P\) was assumed to be known a priori for the simulations with synthetic data. For the simulations with real data, we used the same values as in previous works that studied the same HIs. When such knowledge is not available, this number can be estimated from the HI [89], although care must be taken to avoid overestimating the number of EMs in the scene [90].

To evaluate the results of the algorithms on the synthetic datasets, we considered the normalized mean squared error (NRMSE), defined as \(\text{NRMSE}_{\mathbf{X}}=\|\mathbf{X}-\mathbf{\widetilde{X}}\|_{F}/\|\mathbf{X}\|_{F}\), computed between a matrix \(\mathbf{X}\) and its estimated version \(\mathbf{\widetilde{X}}\), where \(\|\cdot\|_{F}\) denotes the Frobenius norm of a matrix or higher order tensor. We evaluated the NRMSE between the abundances (NRMSE\({}_{\mathbf{A}}\)) and between the EMs (NRMSE\({}_{\mathbf{M}}\)). We also evaluated the normalized image reconstruction error (NRMSE\({}_{\mathbf{Y}}\)); however, we emphasize that a low values of NRMSE\({}_{\mathbf{Y}}\) do not imply good unmixing performance.

We also compute the spectral angle mapper (SAM) between the true and the estimated EMs at each pixel:

\[\text{SAM}_{\mathbf{M}}=\frac{1}{N_{U}}\sum_{n=1}^{N_{U}}\sum_{j=1}^{P}\arccos \left(\frac{\mathbf{m}_{n,j}^{\top}\mathbf{\widetilde{m}}_{n,j}}{\|\mathbf{m}_{n,j}\|_{2} \|\mathbf{\widetilde{m}}_{n,j}\|_{2}}\right). \tag{38}\]

where \(\mathbf{m}_{n,j}\) and \(\mathbf{\widetilde{m}}_{n,j}\) are the true and estimated versions of the \(j\)-th EM in the \(n\)-th pixel, respectively.

### _Synthetic data_

Two synthetic datasets with ground truth were considered, one containing nonlinearity in the mixing process and another containing EM variability. These two nonidealities are intruduced in separate datasets to have a more controlled experimental scenario, particularly since some of the compared methods were devised taking only nonlinearity or EM variability into account. However, in real datasets, especially over large scenes, these effects can occur concurrently.

Data with nonlinearityTo generate the first datacube (DC1), with \(N_{U}=2500\) pixels, we considered synthetic abundance maps displayed in the first row of Figure 3 (left) and \(P=3\) EM signatures with \(L=224\) spectral bands selected from the USGS Spectral Library. The reflectance of each pixel was generated according to the BLMM \(\mathbf{y}_{n}=\mathbf{M}\mathbf{a}_{n}+\sum_{i=1}^{P}\sum_{j=i+1}^{P}a_{n,i}a_{n,j}\mathbf{m }_{i}\odot\mathbf{m}_{j}+\mathbf{e}_{n}\), with \(\mathbf{e}_{n}\) being white Gaussian noise with an SNR of 30dB.

Data with endmember variabilityTo generate the second datacube (DC2), with \(N_{U}=2500\) pixels, we considered synthetic abundance maps displayed in the first row of Figure 3 (right). To incorporate EM variability in this dataset, sets of EM signatures from \(P=5\) pure materials (roof, metal, dirt, tree and asphalt) with realistic variability were first manually extracted from a real HI. Then, to generate each pixel \(\mathbf{y}_{n}\), spectral signatures for each EM were then randomly sampled from these sets and used in a generalized version of the LMM \(\mathbf{y}_{n}=\mathbf{M}_{n}\mathbf{a}_{n}+\mathbf{e}_{n}\) with pixelwise EM matrices, where \(\mathbf{e}_{n}\) was white Gaussian noise with an SNR of 30dB.

Note that the two synthetic datacubes incorporate nonidealities seen in practical scenes. The BLMM, used in DC1, has been experimentally shown to accurately describe multiple scattering in mixtures containing vegetation [106]. The different EM samples used to construct DC2 were manually extracted from a real HI, being representative of the variability encountered in a real image.

DiscussionThe quantitative and visual results for both data cubes are presented in Table I and Figure 3, respectively. It can be seen that ID-Net obtained the best abundance estimation accuracy NRMSE\({}_{\mathbf{A}}\) for both data cubes, with improvements of 39% and 25% over the second-best results. It also obtained a good EM estimation performance, with the best NRMSE\({}_{\mathbf{M}}\) for DC1 and the best SAM\({}_{\mathbf{M}}\) for DC2. The visual results of Figure 3 show that although most methods were able to obtain a rough separation between the different materials in the scene, the reconstructions by ID-Net were the closest to the ground truth. This can be observed most clearly, for instance, in the abundance maps of the Metal EM in DC2. The smallest reconstruction errors NRMSE\({}_{\mathbf{Y}}\) were obtained by NDU, which employs a nonparametric model with many degrees of freedom, which can represent a pixel arbitrarily well. However, this does not imply better unmixing performance (i.e., better abundance and EM reconstructions). The execution times of ID-Net, also shown in Table I, are on the same order of magnitude of the more complex competing algorithms (such as PLMM, DeepGUn and NDU) but considerably higher than the faster ones (such as EGU-Net).
### _Real data_

We also evaluate the performance of the methods on the Samson, on the Jasper Ridge, and on the Cuprite HIs. These images were acquired by the AVIRIS instrument, which samples the spectra on 224 spectral bands. Bands corresponding to water absorption regions or with a low SNR were removed from all datasets, which resulted in 156 bands for the Samson HI, 198 bands for the Jasper Ridge HI, and 188 bands for the Cuprite HI. Previous works have shown that the Samson and Jasper Ridge HIs contain three and four EMs, respectively [30, 32], while fourteen EMs are used for the Cuprite HI [104, 107]. These datasets are relevant since they represent nonidealities such as EM variability, which can appear due to non-trivial variations in illumination from irregular terrain topography and due to the intrinsic variations of the materials such as soil, vegetation and road, and nonlinear mixtures, which are often present in mixtures containing vegetation or water, or in mixtures of minerals. A representation of the real datasets can be seen in Figure 4.

The abundances maps reconstructed by the different algorithms and datasets are presented in Figure 5. For the Samson and Jasper Ridge images, it can be noticed that the methods based on deep-learning (particularly, EGU-Net, RBF-AEC and ID-Net) generally provided a better separation between the different materials in the scene when compared to the remaining algorithms. This can be observed by the smaller confusion in the results by those methods between Water and Soil materials for the Samson HI, and between Road and Soil EMs for the Jasper Ridge HI. For the Cuprite HI (for which only the Musocvite, Buddingtonite, Kaolinite and Chalcedony EMs are shown), a similar behavior is observed for the learning-based methods. The abundances obtained by most of the algorithms generally agree with previous studies on this HI [104, 107], however, for some EMs (the Buddingtonite and Chalcedony) the abundances retrieved by

Figure 4: Visual representation of the real HIs used in the experiments.

Figure 3: True and reconstructed abundance maps for the experiments with synthetic data for data cubes DC1 (left) and DC2 (right).
EGU-Net and RBF-AEC show a different spatial distribution compared to the remaining algorithms. However, pixel-by-pixel comparisons on this scene are hard due to the lack of ground truth [108]. When compared to the competing methods, ID-Net was able to obtain a better abundance reconstruction in regions with more heavily mixed EMs. This can be seen more clearly in the Vegetation and Soil materials in the Samson HI, for which the results of EGU-Net and RBF-AEC identify most pixels as being completely pure. Samples of the EMs estimated by ID-Net are shown in Figure 6. It can be seen that the amount of variability was different for each material and image, being higher in the Samson HI, or for vegetation, compared to the Cuprite HI, or water. It is also instructive to analyze the nonlinearity degree (i.e., the contribution of the nonlinear part of the NN) in the estimated abundance coefficients in (13) for each pixel, defined here as \(\eta_{d}=\|\hat{\gamma}_{\phi}^{a,\mathrm{nlin}}(\mathbf{y},\mathbf{M})\|/\big{(}\| \hat{\gamma}_{\phi}^{a,\mathrm{lin}}(\mathbf{y},\mathbf{M})\|+\|\hat{\gamma}_{\phi}^{a, \mathrm{nlin}}(\mathbf{y},\mathbf{M})\|\big{)}\), which is shown in Figure 7. It can be seen that the overall level of nonlinearity degree varies across the HI. As expected, \(\eta_{d}\) is generally higher in areas with more mixed pixels, and lower in homogeneous areas with predominantly pure pixels. This is particularly seen in areas containing mixtures of vegetation and water in the Samson HI. However, we emphasize that nonlinearity in (13) and EM variability jointly influence the HU results, and separating the effect of each one is not trivial. The quantitative results of the algorithms are shown in Table II. It can be seen that the GLMM and NDU methods obtained the smallest reconstruction errors (NRMSE\({}_{\mathbf{Y}}\)). This occurs since these methods have a high amount of degrees of freedom and can represent HI pixels very closely. However, we note that achieving low values of NRMSE\({}_{\mathbf{Y}}\) do not imply good abundance estimation performance. The execution times of ID-Net were generally higher but still comparable to those by DeepGUn and NDU (1.33 and 1.26 times on average, respectively). This is a limitation of the proposed method, and in future work we will investigate more efficient inference techniques to address it.

The algorithms based on deep learning led to distinct results due to their different frameworks. RBF-AEC and EGU-Net, which are based on AEC NN architectures, obtained abundances that were generally sparse, with many pure components. DeepGUn, which is based on a matrix factorization framework with generative EM models, led to abundances comparatively closer to model-based methods such as the GLMM. ID-Net, which leverages a statistical model of the unmixing processes, gave results with pure components in homogeneous regions and mixed pixels in the regions of transitions between materials. Thus, while the differences between the results of these methods are not straightforward to explain, we can observe how they change as we transition from model-based towards fully neural network-based algorithms.

## VI Conclusions

This paper proposed a hyperspectral unmixing method considering both nonlinearity and EM variability based on a deep disentangled variational inference framework. Both the model parameters and a tractable approximation of the posterior distribution were learned end-to-end by maximizing a lower bound to the likelihoods of supervised and unsupervised data. A self-supervised learning strategy was considered to leverage the benefits of semi-supervised learning algorithms while using only data extracted from the observed image. The abundances and EMs were disentangled through the use of appropriate independence assumptions on the variational posterior. Moreover, an interpretable model was designed by using unrolled optimization-based and two-stream (linear/nonlinear)

Figure 5: Reconstructed abundance maps for the experiments with real data for the Samson (left), Jasper Ridge (middle) and Cuprite (right) HIs.
NN architectures, allowing for an adjustable amount of nonlinearity in the model. Experimental results on both synthetic and real datasets showed that the proposed method can achieve state-of-the-art performance with execution times comparable to that of algorithms such as NDU and DeepGUn. Several interesting directions for future work can be considered, such as developing more efficient inference strategies that lead to algorithms with smaller computational cost; designing self-supervised learning approaches to extract the training data from the HI which do not depend on the presence of multiple pure pixels; and exploiting the statistical dependence among abundances of spatially adjacent pixels.
