Collaborative sparse regression using spatially correlated supports - Application to hyperspectral unmixing

This paper presents a new Bayesian collaborative sparse regression method for linear unmixing of hyperspectral images. Our contribution is twofold; first, we propose a new Bayesian model for structured sparse regression in which the supports of the sparse abundance vectors are a priori spatially correlated across pixels (i.e., materials are spatially organized rather than randomly distributed at a pixel level). This prior information is encoded in the model through a truncated multivariate Ising Markov random field, which also takes into consideration the facts that pixels cannot be empty (i.e., there is at least one material present in each pixel), and that different materials may exhibit different degrees of spatial regularity. Second, we propose an advanced Markov chain Monte Carlo algorithm to estimate the posterior probabilities that materials are present or absent in each pixel, and, conditionally to the maximum marginal a posteriori configuration of the support, compute the minimum mean squared error estimates of the abundance vectors. A remarkable property of this algorithm is that it self-adjusts the values of the parameters of the Markov random field, thus relieving practitioners from setting regularization parameters by cross-validation. The performance of the proposed methodology is finally demonstrated through a series of experiments with synthetic and real data and comparisons with other algorithms from the literature.


I. INTRODUCTION
Spectral unmixing (SU) of hyperspectral images is a challenging problem that has received a lot of attention over the last few years [1]- [3]. It consists in identifying the materials (endmembers) present in an image and simultaneously quantifying their fractions or proportions within each pixel (abundances). This source separation problem has been widely studied for applications where pixel reflectances are linear combinations of pure component spectra [4]- [8]. The typical SU processing pipeline is then decomposed into three main estimation steps: the estimation of the number of different materials present the image, the estimation or extraction of their spectral signatures, and finally the quantification of their abundances.
Abundance estimation is known to be a challenging problem, in particular in scenarios involving materials with similar spectral signatures. In these cases, exploiting prior knowledge about the problem can improve estimation performance dramatically. A particularly important form of prior knowledge is that the number of materials within each pixel is typically much smaller than the total number of materials present in the scene (this property is accentuated in modern images that are acquired with high spatial resolution sensors). In other words, the abundance vectors are generally sparse. Sparsity also arises naturally when SU is performed with a dictionary or library containing the spectral signatures of a large number materials, some of which are possibly present in the scene [9].
Once sparsity is taken into consideration, SU can be conveniently formulated as a sparse regression problem whose objective is to jointly identify the materials within each pixel and to quantify their abundance. This regression problem is often solved by penalised maximum likelihood estimation, which can be efficiently computed with state-of-the-art optimisation algorithms (typically an 1 penalty is used to promote sparse solutions) [10]. Recently, Iordache et al. [11] proposed a collaborative sparse regression technique (CLSunSAL) based on an 2,1 penalty function that enforces group-sparsity for the abundances of each material.
This method was further improved in [12] by introducing a pre-processing step that identifies the elements from the spectral library that are more likely present in the image. The resulting MUSIC-CSR algorithm solves a sparse regression problem that is collaborative in the sense that all the image pixels are used to identify the active endmembers. Sparse regression for SU can also be successfully performed within the Bayesian framework. For example, Dobigeon et al. [13] and Eches et al. [14] propose Bayesian models and Monte Carlo algorithms to identify the active endmembers in an HSI from a spectral library while ensuring that the abundances of absent endmembers are zero. Note that library-based methods are not the only strategy to address the absence of pure pixels (see [4], [15]- [17] for more details).
It is widely acknowledged that collaborative sparse regression methods can produce very accurate SU results. Collaboration is key because it improves the estimation of the support of the sparse abundance vectors, thus reducing significantly the number of unknowns. However, most existing collaborative techniques only exploit global information and therefore can only seek to determine if an endmember is present/absent in the entire image (that is, can only estimate the union of the supports of all the abundance vectors in the image). As a result, global collaborative techniques may overestimate significantly the support of the actual abundance vector of each pixel. This paper presents a new collaborative sparse regression technique that exploits the local spatial correlations in the image to accurately detect the endmembers that are active/inactive in each pixel. Precisely, we present a Bayesian model that simultaneously promotes sparsity on the abundance vectors, and spatial correlation between the abundance vectors' supports (i.e., non-zero elements), modelling the spatial presence and absence patterns of materials in the scene. This approach differs significantly from the strategies adopted in the previous works [18]- [22], where spatial correlations are introduced by regularising the abundance values or the nonlinear effects occurring in the image. The latter strategies perform well in hyperspectral images with low spatial resolution, and composed mainly of homogenous regions, but are inadequate for high resolution images and for images involving complex scenes, small targets, and textures or fluctuations in the abundances. By operating directly at the level of the abundance vectors' support, our model is able to capture spatial correlations in a more subtle manner and produce accurate estimation results in challenging scenes involving, for example, crops (textured scenes) and isolated trees.
where a r,n is the mixing coefficient associated with the rth endmember in the nth pixel and a n = [a 1,n , . . . , a R,n ] T , Eq. (1) can be conveniently expressed in matrix notation as y n = Ma n +e n . This paper considers the supervised spectral unmixing problem of the hyperspectral image Y, i.e., the estimation of the R × N abundance matrix A = [a 1 , . . . , a N ]. More precisely, we consider Bayesian methods for estimating A given Y (and the endmember matrix M) subject to the following two sets of physical constraints: first, the abundances are non-negative quantities, i.e., a r,n ≥ 0 ∀r, n; second, there is at least one material present in each pixel and therefore a n 0 > 0, where · 0 denotes the 0 vector pseudo-norm. We also assume that the values of the noise variances σ 2 are unknown, though prior knowledge, if available, can be easily integrated into the model.
In a manner akin to [23], we model explicitly the sparsity of A by using the decomposition where Z ∈ {0, 1} R×N is a matrix of Bernoulli variables that "labels" each material as present (active) or absent (inactive) in each pixel, X ∈ R R×N is a matrix with positive entries that (jointly with Z) quantifies the abundances, and denotes the Hadamard (term-wise) matrix product. This decomposition is particularly useful in Bayesian sparse regression problems because it allows eliciting separate statistical models for a n 's support (through modelling Z), and for the values of the positive elements of a n (through X).
The next section presents a Bayesian model for estimating Z and X subject to the physical constraints discussed above. A key aspect of this model is that it will capture the fact that the pixels in which a material is active (or inactive) generally form spatial clusters (i.e., exhibit spatial group sparsity). In difficult unmixing scenarios, exploiting this strong prior information can significantly improve estimation results, as will be shown in this paper. lihood and the prior distribution of the parameters of interested Z and X, as well as for the other unknown parameters in the model (e.g. σ 2 ) that will be subsequently removed by marginalisation (i.e., integrated out of the model's joint posterior distribution).

A. Likelihood
From the observation model (1) and the parametrisation of A described in (2), the likelihood of the image Y given the unknown parameters is where z n (resp. x n )is the nth column of Z (resp. X) and p N (y n |M(z n x n ), Σ 0 ) is the probability density function of a multivariate Gaussian vector with mean vector M(z n x n ) and diagonal covariance matrix Σ 0 .

B. Prior distribution of Z
As explained previously, a key aspect of the proposed Bayesian model is to take into account the fact that the pixels in which a given material is present or absent typically form spatial groups or clusters (as opposed to being randomly distributed in space). From a modelling viewpoint, this can be represented by correlating the Bernoulli variables or labels z r,n across the spatial dimension indexed by n. This can be achieved, for example, by stating that if a certain material is present (or absent) in a given pixel, this increases the probability of its presence (or absence) in neighbouring pixels. Taking into account that each material can exhibit its own spatial configuration, and the constraint that there must be at least one material present in each pixel, we propose to assign Z the following truncated multivariate Ising Markov random field prior with β = {β 1 , . . . , β R }, and where δ(·) is the Kronecker delta function and V(n) denotes the set of neighbours of pixel (n) (in this paper we have used the 8-pixel neighbourhood). The hyper-parameters β 1 , . . . , β R act as regularization parameters that control the degree of spatial smoothness or regularity associated with each endmember, accounting for the fact that different materials may exhibit different spatial distributions.
To gain intuition about the proposed prior, we note that setting β = 0 in (5) and considering the prior (9) for the matrix X leads to a Bernoulli-Gaussian type prior for the abundances which is closely related to the 0 -2 penalty often used for sparse regression. In the general scenarios where β = 0, Eq. (5) introduces spatial correlations between the components of Z and, together with (9), leads to an 0 -2 -type penalty promoting spatial group sparsity for the abundances. Prior (5) can also be understood as a collaborative prior, in the sense that, by assigning higher probabilities to configurations in which the supports of the abundance vectors a n are spatially correlated, it pools or shares information between the sparse regressions that take place at each pixel. Crucially, by capturing the spatial correlations that occur naturally in hyperspectral images, this prior can improve significantly estimation results.

C. Prior distribution of X
We assign the elements of X the following hierarchical prior distribution parametrised by some fixed hyper-parameters (γ, ν), and where we note that the prior on x r,n |s 2 r is truncated to R + to reflect the positivity of x r,n . This prior is very flexible and can be adjusted to represent a wide variety of prior beliefs (in all our experiments we used γ = 2.1 and ν = 1.1, corresponding to a weakly informative prior for s 2 r with 80% of its mass in [0, 1].) An important advantage of the hierarchical prior (9) is its natural capacity to encode prior dependences between the abundances. We expect the abundance coefficients associated with the same material to exhibit correlations, in particular in terms of their scale. This belief is encoded in (9) by defining a common parameter s 2 r for each material or endmember m r , which is shared by all the abundances related to that material. Therefore, the hierarchical structure of (9) operates as a global pooling mechanism that shares information across the rows of A (i.e, the abundance coefficients associated to the rth material) to improve estimation accuracy.
Assuming an exchangeable structure where the abundances are prior independent given the hidden variables s 2 = [s 2 1 , . . . , s 2 R ] T , we obtain the following following joint prior for X, s 2 with f (X|s 2 ) = r,n f (x r,n |s 2 r ) and f (s 2 ) = r f (s 2 r |γ, ν). Also notice that by using the hierarchical structure (9) we obtain conjugate priors and hyper-priors for x r,n and s 2 r ; this leads to inference algorithms with significantly better tractability and computational efficiency, which is crucial given the high dimensionality of X.
Finally, notice that in the proposed model the spatial dependences in the abundance maps are encoded at the level of the abundance supports (though the prior on Z), and not directly through the values of the abundances A = Z X as it is the case in some previous models (see for example [19]). The motivation for this modelling choice is that modern high-resolution hyperspectral images often exhibit textures and fine detail that are not well described by models that promote smooth or piecewise-constant abundances maps. The experiments reported in this work show that the two approaches to modelling spatial correlation have complementary strengths and weaknesses, suggesting that to further improve estimation results future models should consider both levels of spatial dependences simultaneously.

D. Prior distribution of the noise variances σ 2
In this paper we consider that there is no significant prior knowledge available about the values of the noise variances and assign each σ 2 its non-informative Jeffreys prior [24] f (σ 2 ) ∝ 1 Note that in scenarios where prior knowledge about σ 2 is available, this can be easily integrated into the model by replacing (12) with an inverse gamma conjugate priors with hyper-parameter values reflecting this prior knowledge.

E. Regularisation parameter β
A main advantage of Bayesian methods is that they allow estimating the appropriate amount of regularisation from data, thus freeing practitioners from the difficulty of setting regularisation parameters by cross-validation. Indeed, there are several Bayesian strategies for selecting the value of the regularisation parameter β in a fully automatic manner (see [25] for a recent detailed survey on this topic). In this paper we use the empirical Bayes technique recently proposed in [26], where the value of β is estimated by maximum marginal likelihood. Using Bayes' theorem, and taking into account the conditional independences of the Bayesian model (see Fig. 1), the joint posterior distribution of Z, X, σ 2 and s 2 given the observations Y, the library of spectral signatures M and the model's fixed parameters β, γ and ν, is given by The following section presents a Monte Carlo algorithm to perform Bayesian inference using the proposed model.

IV. BAYESIAN INFERENCE USING GIBBS SAMPLING
The Bayesian model defined in Section III specifies a joint posterior distribution for the unknown parameters Z, X, σ 2 , s 2 given the fixed quantities Y, M, γ, ν and the hyperparameter β which is unknown but represented as a deterministic parameter (whose value will be tuned during the inference procedure). According to the Bayesian paradigm, this posterior distribution fully describes the information about the unknowns that is provided by the data and by the prior knowledge available. However, for unmixing applications it is necessary to summarise this posterior distribution in the form of point estimates; that is, to assign specific values for the unknown quantities of interest (in our problem the abundance vectors). Here we consider the following coupled Bayesian estimators that are particularly suitable for sparse regression problems: the marginal maximum a posteriori (MMAP) estimator [27], [28] for the support of the abundance vectors or "presence maps" and, conditionally on the estimated supports, the minimum mean square error estimator of the abundances where f (z r,n |Y, M, β, γ, ν) = f (Z, X, σ 2 , s 2 |Y, M, β, γ, ν)dZ \zr,n dX, dσ 2 ds, with the matrix Z \zr,n containing the remaining elements of Z once z r,n has been removed and where E [·] denotes the expectation with respect to the conditional marginal density f (x r,n |z r,n , Y, M, β, γ, ν) = f (Z, X, σ 2 , s 2 |Y, M, β, γ, ν)dZ \zr,n dX \zr,n , dσ 2 ds f (z r,n |Y, M, β, γ, ν) .
Note that the abundance estimator (15) is sparse by construction (i.e., E [x r,n |z r,n = 0, Y, M, β, γ, ν] = 0) and that, by marginalising out the other unknowns, it automatically takes into account the uncertainty about σ 2 and s 2 .
The choice of specific Bayesian estimators to summarise the posterior distribution is a decision-theoretic problem that depends on the model and the application considered [27]. In the model described in this work, the two quantities of interest Z and X are very different in nature, and as a result their posterior information is best summarised with different estimators.
The labels Z are binary variables that describe a quantitative aspect of the model; that is, they parametrise the hypotheses that materials are present or absent in each pixel. Selecting the hypotheses with highest posterior probability leads to marginal MAP estimation [27], which in this case operates as a model-selection tool. Moreover, the conditional posterior distribution of X|Z represents the uncertainty regarding the abundance values for a specific model or configuration (in particular the one identified by marginal MAP estimation). To summarise this posterior distribution we use the MMSE estimator because it is optimal with respect to any quadratic loss function, and approximately optimal with respect to any convex loss function [27]. Furthermore, in cases where there is strong prior knowledge justifying the use of a sum-to-one constraint on the abundances, this information can be incorporated to the inferences by using an MMSE estimator constrained to the simplex. However, we have observed that this property does not hold on our images, in part because of the effects of mild shadows in the scene.
Computing (14) and (15) is challenging because it requires having access to the univariate marginal densities of z r,n and the joint marginal densities of (x r,n , z r,n ), which in turn require computing the posterior (13) and integrating it over a very high-dimensional space. Fortunately these estimators can be efficiently approximated with arbitrarily large accuracy by Monte Carlo integration. Precisely, it is possible to compute (14) and (15) by first using a Markov Chain Monte Carlo (MCMC) computational method to generate samples asymptotically distributed according to (13), and subsequently using these samples to approximate the required marginal probabilities and expectations.
The remainder of this sections provides details about the main steps of the proposed Gibbs sampler, termed Collaborative sparse Unmixing (CSU) and summarised in Algo. 1 below.

A. Sampling the label matrix Z
The label matrix Z is updated pixel-wise by iteratively simulating from the distribution of the labels at each pixel given the other pixels. Precisely, the distribution of the vector z n given the matrix Z \zn containing the remaining elements of Z once z n is removed is given by whereỹ n = y n − M(x n z n ), for ||z n || 0 > 0 and P (z n |Y, Z \zn , X, σ 2 , s 2 ) = 0 otherwise.
Algorithmically, this simulation step can be achieved by indexing the 2 R − 1 admissible configurations of z n and then randomly selecting a specific one with probability defined in (18).

B. Sampling X
The conditional distribution of X given the other unknown parameters can be factorised pixel-wise as a product of N marginal distributions that can be efficiently sampled independently and in parallel where Ω = (R + ) R is the positive orthant of R R and and where S = diag(s 2 ) and D n = diag(z n ) are diagonal matrices with diagonal elements given by s 2 and z n . For completeness, the derivation of 21 is provided in Appendix. In this paper we use the method [30] to simulate efficiently from (21).

C. Sampling the noise variances σ 2
It can be easily shown that the noise variances are (conditioned on the other parameters) a posteriori dependent and can thus be updated in a parallel manner. Precisely, the conditional distribution associated with σ 2 has a simple closed form expression and is given by with E σ = ỹ ,: that can be easily sampled independently and in a parallel manner

V. VALIDATION WITH SYNTHETIC DATA
This section demonstrates the proposed methodology on a series of experiments conducted using synthetic data. An applications to a real hyperspectral image is reported in Section VI.

B. Supervised unmixing
In this first experiment we consider that the materials present in the images I 1 and I 2 are perfectly known and we estimate the abundance vectors with CSU, CLSunSAL and NCLS. The CSU algorithm has been implemented using N MC = 3000 and N bi = 1000, and by allowing the algorithm to self-adjust the regularisation parameter β with the technique proposed in [26]. of-the-art sparse regression algorithm SunSAL [10] and the SunSAL algorithm using the total variation regularization (SunSAL-TV) [19], whose regularisation parameter values are adjusted to provide the best abundance estimates (in terms of RMSE) for each scenario.
For completeness we also compare with the widely used Non-negatively Constrained Least-Squares algorithm (NCLS) [5], which solves a maximum likelihood problem and does not exploit any prior information about the abundance vectors.
The associated computation times for a Matlab implementation on a 3.0GHz Intel Xeon quad-core workstation are provided in Table I. The estimated supports obtained with CSU, SunSAL, SunSAL-TV and NCLS for I 2 (lowest SNR) are depicted in Fig. 3. The presence maps associated with I 1 are similar and are not presented here due to space constraints. For SunSAL, SunSAL-TV and NCLS, the detection maps have been obtained by thresholding the estimated abundances with a threshold arbitrarily set to ρ = 0.01. We observe that the results obtained with CSU are in good agreement with the ground truths. On the other hand, the abundances obtained with SunSAL, SunSAL-TV and NCLS are significantly less accurate, thus confirming that taking into account the spatially correlations between the supports abundance vectors is key to achieving accurate estimation results in this scenario.
This valuable prior knowledge is all the more important in low SNR scenarios such as the one depicted in Fig. 3.
For numerical comparison, we computed the root mean square error (RMSE) RMSE n = â n − a n 2 , that quantifies the average accuracy of the estimated abundancesâ n with respect to the truth a n at the n-th pixel. We also consider the abundance angle distance (AAD) AAD n = cos −1 â T n a n â n 2 a n 2 , which is not sensitive to scaling factors between actual and estimated abundance vectors. The lower the RMSEs and AADs, the better the abundance estimation performance.
The first five rows of Table II show the average RMSEs and AADs, for CSU, SunSAL, SunSAL-TV, NCLS, and the oracle NCLS (o-NCLS), which consists of applying NCLS using only the active materials in each pixel (i.e., with perfect knowledge of the support of the abundance vectors). We observe that CSU provides significantly more accurate estimations than SunSAL, SunSAL-TV and NCLS, achieving average RMSEs and AADs that are close to the oracle. Again, the good performance of CSU results from the collaboration between pixels introduced by the prior model (5). By comparing the first and second columns of Table   II we confirm that this prior knowledge becomes all the more important as the noise level increases. These results show that using only sparsity (SunSAL) does not lead to significant improvements, perhaps because the 1 regularisation is not appropriate for our images, or because it would require adapting the regularisation parameters to each endmember (this is done automatically in CSU). We also observe that SunSAL-TV does not achieve significantly better results, perhaps due to the regions where the abundances exhibit fluctuations that are not well modelled by the TV. As mentioned previously, when both the abundance supports and the abundance values exhibit significant spatial dependencies, one should consider methods that account for this prior information (e.g., the SunSAL-TV method if the abundance maps are expected to be piecewise constant across the image). However, in many highresolution images and regions in images, only the abundance supports exhibit significant spatial dependencies; in this case regularising the abundance values may lead to over-smooth estimation results. In this paper we show that by modelling these dependencies it is possible to improve unmixing performance with respect to methods that assume that abundances are uncorrelated.
Finally, in order to further highlight the good performance of the proposed prior model we have also computed the reconstruction error (RE) defined as RE n = Mâ n − y n 2 .
Note that this error essentially measures the likelihood of the abundance estimates given the observed image Y and assuming the noise is identically distributed in the L spectral bands.
However, because the SU problem is not well-posed, the capacity of the likelihood to identify good solutions is severely limited and the additional information provided by the prior model is key to deliver accurate estimation results. In these scenarios Bayesian methods, which combine observed and prior information, can greatly outperform other estimation techniques.
The average RE for each method is reported in Table III. We observe that all the algorithms considered exhibit very similar REs. However, we know from Table II that CSU outperforms significantly the other methods in terms of abundance estimation accuracy. The contrast between these two figures of merit confirms that the superior performance of CSU is directly  related to the proposed prior model, which captures the correlations between the supports of the abundance vectors and effectively introduces a means of collaboration that allows sharing information between pixels and increasing robustness to noise.

C. Semi-supervised unmixing
We now consider that two additional endmembers are incorrectly included in the library M, although they are not present in the scene. To make the SU problem particularly challenging, the two additional endmembers are Olivine 2 and Adularia, whose spectral signatures are  Finally, the five bottom rows of Table II show the RMSEs and AADs obtained with CSU, SunSAL, SunSAL-TV, CLSunSAL and NCLS for I 1 and I 2 with R = 7. We observe that

VI. APPLICATION TO REAL A HYPERSPECTRAL IMAGE
This section presents an application of the proposed CSU method to a real hyperspectral image acquired by the Hyspex hyperspectral scanner over Villelongue, France (00 • 03 W and 42 • 57 N). This images was acquired in 2010 as part of the Madonna project, and is composed of L = 160 spectral bands covering from the visible to near infrared spectrum and with a spatial resolution of 0.5m (for more details about the data acquisition and pre-processing steps see [33].). This dataset has previously been studied in [22], [33], [34] and is mainly composed of forested and urban areas. Due to the high spatial resolution of the image, materials present in a pixel are likely to be also present in its neighbour pixels; that is, we expect significant spatial correlations between the mixture supports. Here we evaluate the proposed unmixing method on the region of interest of size 100 × 100 pixels depicted in Fig. 6. This region is composed mainly of trees and grass, with R = 6 endmembers related to soil, two types of grass, two types of trees, and an additional endmember modelling attenuation effects mainly related to shade. The spectral signatures for these endmembers have been extracted manually from the data by using our prior knowledge about the scene [33]. Note that for this image and within the LMM framework, the attenuation effects have to be modelled as an additional endmember because they operate differently on the different spectral bands of the reflectance spectra (otherwise, if the attenuations acted as a scaling factor affecting all the bands similarly, we would be able to handle this by relaxing the abundance sum-to-one constraint). Fig. 7 shows the presence maps for each material estimated with CSU (using N MC = 5000, N bi = 1000 and by allowing the algorithm to self-adjust the regularisation parameter β with the technique [26]) and with NCLS, SunSAL and SunSAL-TV (whose parameters have been optimised in a fully supervised manner to obtain the sparser and/or smoother abundance maps without degrading significantly the reconstruction error, and by using a detection threshold of ρ = 0.01). Fig. 8 shows the total number of materials present in each pixels, as estimated by each method. Lastly, Fig. 10 shows the abundance maps estimated with each method. These figures clearly show that CSU provides sparser and spatially smoother presence (support) maps than the other methods. We also observe that, as expected, SunSAL-TV provides smoother abundance maps than the other methods (see Fig. 10). This property of SunSAL-TV is beneficial in regions where the materials are absent; however, where materials are present, it may smooth out fine detail (e.g., textures). Since no abundance ground truth is available for this data set, it is difficult to determine which algorithm provides the more accurate results. However, due to the high spatial resolution of the image and the structure of the materials considered, it is reasonable to assume that abundances will exhibit some degree of local heterogeneity, and that CSU and SunSAL-TV will outperform each other in different ways. Moreover, from Figs. 7 and Fig. 10 we also note that CSU detects a higher lever of attenuation effects than SunSAL and SunSAL-TV, particularly in pixels containing vegetation where strong shadowing effects occur.
Finally, as a model checking analysis, we plot in Fig. 9 the marginal noise variances estimated by CSU (i.e., diagonal elements of the full estimated covariance matrix). For comparison we also include the estimates obtained with Hysime [35]. We observe that the estimation results produced by CSU are very smooth and physically realistic, suggesting good fit to data (recall that CSU assumes that the noise variances are prior independent, hence the smoothness of these posterior estimates arises from the data and from correlations with other model parameters). We also note that, although the estimation results produced by CSU and Hysime are different, the noise variances follow similar profiles. Lastly, for completeness Table IV reports the reconstruction errors associated with CSU, NCLS, SunSAL and SunSAL-TV. We observe that all methods achieve comparable reconstruction errors, with NCLS and SunSAL having a slightly lower reconstruction error than CSU and SunSAL-TV.
This difference is likely due to the fact that NCLS and SunSAL do not enforce spatial regularity, and are hence more prone to overfitting the data.  Due to its computational complexity, the algorithm presented in this paper can only be directly applied to problems with small numbers of endmembers (e.g., R ≤ 25). For problems with larger libraries it is computationally more efficient to generate samples from a relaxed posterior in which the "non-empty pixel" constraint (7) of the MRF is removed, and then reintroduce this constraint by importance sampling. Another possibility for problems with large libraries is to use the MUSIC-CSR algorithm [12] as a pre-processing step to identify endmembers that are absent from the scene, and then apply our method using a pruned library.
As mentioned previously, previous works have considered that the values of the abundance vectors are spatially correlated, whereas we considered correlations between their supports (i.e., material presence and absence patterns). Investigating new models that exploit both types of correlations is an important perspective for future work. For many applications or hyperspectral datasets, it makes sense to consider additional abundance constraints, such as the sum-to-one constraint to further improve unmixing results. Embedding this constraint within 0 -type sparse regression models is a challenging problem that is currently under investigation. Another perspective for future work is to investigate more sophisticated spatial models that describe hyperspectral images more accurately, in particular in complex scenes with numerous materials and non-linear effects. Finally, we believe that the methodology presented in this paper could be interesting for other regression problems that exhibit structured sparsity and intend to further investigate this is the future.

ACKNOWLEDGMENTS
The authors would like to thank Prof. Jean-Yves Tourneret and Dr Nicolas Dobigeon, from the University of Toulouse, IRIT-ENSEEIHT, France, for interesting discussion regarding this work.

APPENDIX: ON THE CONDITIONAL DISTRIBUTION OF THE ABUNDANCE VALUES
Due to the conjugacy of the priors (9) and the likelihood (3), the posterior distribution f (x n |y n , z n , σ 2 , s 2 ) ∝ f (y n |x n , z n , σ 2 )f (x n |s 2 ) is a multivariate Gaussian distribution restricted to Ω, i.e., the positive orthant of R R . Moreover,∀x n ∈ Ω log (f (x n |y n , z n , σ 2 , s 2 )) = log (f (y n |x n , z n , σ 2 )) log (f (x n |s 2 )) , where S = diag(s 2 ) and D n = diag(z n ) are diagonal matrices with diagonal elements given by s 2 and z n and c 1 , c 2 are real constants (independent of x n ). By identification, we obtain