IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. ?, NO. ?, 2 II. OPTIMIZATION STRATEGIES A. Technical Optimization (TO) is introduced. Its influence on the separation performances is also discussed. The issue of estimating the number of sources, or “intrinsic dimension” [28], will not be addressed in this article. Several methods have been proposed in the literature [29], [30]. This paper is organized as follows. Section II describes the proposed implementation strategies adopted for this work. Sec- tion III summarizes the improvements related to the technical aspects of memory storage and computation issues. Section IV discusses the performances of the resulting algorithms when the pixel selection preprocessing step is introduced. strained Endmembers) has been proposed in [12] to minimize the size of the simplex formed by the estimated endmembers. However, as noted in [13], results provided by ICE strongly depend on the choice of the algorithm parameters. More re- cently, Jia and Qian have developed in [14] complexity-based BSS algorithms that exploit the pixel correlations to recover endmember signatures. In [15], Miao et al. have introduced a non-negative matrix factorization (NMF) algorithm with an additivity penalty on the abundance coefficients. Similarly, other constrained NMF approaches exploiting smoothness and sparseness features have been considered in [16]. Note that all the strategies described above are based on an optimization scheme to minimize a penalty criterion. Consequently, they may suffer from convergence issues, e.g., due to the presence of local maxima and the large number of parameters to be estimated. The optimization consists of two independent parts which will be referred to: (i) Technical Optimization (TO) to reduce the memory footprint, the average cost of algorithmic oper- ations, and make smarter reuse of memory (ii) Convex Hull Optimization (CHO) to reduce the number of spectra to be processed. Both parts enabled us to analyze hyperspectral images that so far were not open to analysis. The authors stress that the techniques applied in (i) do not alter the results of the original algorithm (see section III). On the other hand, the optimization strategy (ii) only selects a subset of the original input and therefore may change the results. Impact of the strategy (ii) needs to be evaluated, which will be presented in section IV. Alternatively, the constrained separation problem can be conveniently addressed in a Bayesian framework. Two algo- rithms that perform unsupervised separation under positivity and sum-to-one constraints have been recently proposed [17], [18]. These algorithms are based on hierarchical Bayesian modeling to encode prior information regarding the obser- vation process, the parameters of interest and include the positivity and full additivity constraints. The complexity of the inference from the posterior distribution of the parameters of interest is tackled using Markov chain Monte Carlo (MCMC) methods [19], [20], which has been proposed to analyze hyperspectral images [21]. The algorithmic details are not described here. The reader is invited to read the pseudo- codes summarized in Algo.’s 1 and 2, and to consult [17] and [18] for further details. The only difference between the two BPSS algorithms is the sampling of the abundance vectors a$_{p}$ ( p = 1 , . . . , P ). However, since these algorithms rely on MCMC methods, the computation time drastically increases with the image size and these algorithms have not been applied for large scale data processing in spite of their high effectiveness. The algorithms introduced in [17], [18] and referred to as BPSS and BPSS2, respectively, could be successfully launched on an image of a restricted size, typically of a few thousand pixels. The main goal of this work is to optimize the memory requirement of these algorithms to process a whole hyperspectral image of 100000 spectra as it typically occurs in Earth and Planetary Science. Since the time requirements of the computation increase drastically for a larger number of pixels and a larger number of sources, another challenging objective is to reduce as much as possible the computation time. In that respect, our proposal is to discuss the memory storage, the data representation, the operating system archi- tecture and the computing parallelization. These algorithms have been implemented in MATLAB c ○ for this work but future implementations will be done in other languages as well. 1) Memory: Thanks to the MATLAB c ○ profiler, it can be noticed that the main limitation of the BPSS implementation is the contiguous memory. Fragmentation may occur when variables are resized after memory allocation. In this case, the memory management might not be able to allocate a chunk of memory that is large enough to hold the new variable. Significant garbage collection may set in, which may have a significant performance impact. In our case, to reduce the impact of garbage collection, pre-allocating the matrices and work with global variables has been found to be useful. The aim of this article is to discuss some implementation strategies which allow one to apply these algorithms to real hyperspectral data, even if images are large. Previous works about blind source separation of hyperspectral images have been proposed [22]–[24] but only few use positivity/sum- to-unity constraints [25]. To overcome this difficulty, a first approach has been proposed in [25] to combine Independent Component Analysis (ICA) and Bayesian Positive Source Separation (BPSS). Firstly, applying an ICA algorithm (such as JADE [26] or FastICA[27]) to the hyperspectral images is applied to get a rough spatial classification of the scene and to sample relevant pixels (i.e., from each class , the pixels whose spectra are mostly uncorrelated are selected). Secondly, the spectra associated to these pixels will serve in the Bayesian separation algorithm to estimate the endmember spectra. Fi- nally, the abundances can then be estimated on the whole image using the estimated spectra. However, this strategy presents a limitation related to the difficulty to determine the number of pixels to retain from each independent component class. In this paper, another pixel selection strategy based on the computation of the convex hull of the hyperspectral data 2) Precision: MATLAB c ○ by default computes on double precision. However computing with single data type saves a lot of computation time while providing sufficient arithmetic precision. It has been estimated to win up to 60% computation