Soazig Corbel
email: s.corbel@brgm.fr

Joël Kaiser

Serge Vicentin

Coal mine flooding in the Lorraine-Saar basin: experience from the French mines

Keywords: hard coal, mine flooding, water quality, Lorraine-Saar basin

Coal mining in the Lorraine-Saar basin started in the nineteenth century. The coal deposit was mined both in France and in Germany. On the French side, mining ceased in 2004 after 150 years of activity, and as a consequence, the French coal mines got progressively flooded. This study presents the results of the entire monitoring plan for the different mine water reservoirs of the Lorraine coal basin. The actual speed of mine flooding is compared to prediction studies, mine water quality is detailed and content evolution of iron, manganese and suspended solids is discussed for the different mine reservoirs.

Introduction

Coal mining in the Lorraine-Saar basin started in the nineteenth century. The last French coal was mined in April 2004, and once mining stopped, the water that was infiltrating from groundwater to underground mine workings was no longer pumped out to the surface by the mine operator, Charbonnages de France (CdF). From 2006, mine workings got progressively flooded, creating mine water reservoirs (Figure 1). Before its dissolution in 2007, CdF planned and financed remediation operations, including studies to understand mine flooding and its consequences. After the transfer of responsibility for CdF's facilities to the French State in 2008, the Mine Safety and Risk Prevention Department of the French Geological Survey has been managing this remediation activity.

Mine flooding

In 2001, the mine operator CdF started investigating the mechanisms of mine flooding and its consequences. Although coal mines were still operating in the Lorraine basin, planning ahead for post-mining remediation turned out to be necessary.

Mine flooding predictions

Mine flooding assessment was carried out in two phases: volume computation of the residual mining voids and simulation of mine flooding. The first phase started in 2002 and ended in 2005, with a team dedicated only to that task. First, all mine workings maps were digitized, including coal seams and galleries of all exploitation levels for the different exploitation fields. In 2002, CdF used an isolated exploitation field as a case study, to determine coefficient for residual voids depending on the depth of the mine workings and the backfilling methods (Degas 2002). Then from 2003 to 2005, residual mining voids were computed using the residual voids coefficients and the digitized maps, block after block, level after level, for all the mine workings.

The simulation of mine flooding was done in 2006. At that time, CdF was still pumping out water from the underground mine workings. The flow rates of each mine block were assessed and used for the simulation. They were considered equivalent to the flow rates of infiltrating water from the above groundwater. CdF also assessed the state of all galleries linking exploitation fields to each other, concluding that the main galleries were still open. At last, using the residual mining voids, the flow rates of infiltrating water of each block and the connections between each exploitation field, CdF used a simple computation to make mine flooding predictions.

Predictions versus reality

Mine flooding predictions were made for each exploitation field. The example of the Forbach exploitation field, in the Central-Eastern mine water reservoir, is shown on Figure 2. Two mine flooding scenarios were studied by CdF. The first scenario considered that infiltrations would stay constant until all the mine workings would be flooded, up to the semi-permeable geological layer of the Permian. The second scenario considered that the infiltrations would decrease once the mine flooding reached the last hundred meters of the mine workings.

3

The mine flooding predictions made all over the Lorraine basin turned out to be pretty accurate. For the Forbach exploitation field, the speed of flooding for the first 500 meters was overestimated by only a month (Figure 2), probably due to a slight inaccuracy of the residual voids coefficient. Out of the two scenarios studied by CdF, the scenario taking into account a decrease in the infiltrations is the most realistic. However, the speed of flooding for the last 150 meters before reaching the Permian layer was also overestimated, by just over a year. As the shallow mine workings date from early 1900's, the residual mining voids were not easy to assess. The inaccuracy of the prediction of the last 150 meters is most probably due to an inaccurate residual voids computation.

Feedback on mine flooding predictions

Although CdF used a straightforward way to predict mine flooding for the coal mines of the Lorraine basin, the predictions were quite accurate. The key of mine flooding predictions turned out to be data gathering, and not sophisticated computations. Only staff members of the mine operator were able to properly assess residual mining voids and the state of the galleries connecting exploitation fields. However, even with the experience, those miners had trouble assessing old mine workings, and as discussed previously, this inaccuracy immediately lead to an overestimation of the speed of mine flooding.

Mine water quality

To prevent mine water from contaminating the above groundwater once the mine water level has reached the Permian layer, three mine water pumping stations associated to treatment schemes have been built on the Lorraine coal basin. Those pumping stations enable to keep mine flooding under control. They are run by the BRGM as part of their post-mining missions.

Mine water quality monitoring

Water sampling and water analysis are run on the three pumping stations. The pumping stations are installed in old mine shafts, with access to the mine workings. When the coal mines closed, some of the underground mining facilities could not be removed. After being emptied and cleaned, some of the facilities stayed in place, including machinery, electrical and hydraulic installations. Despite all the precautions taken, some industrial toxic products such as oils and resins might have remained underground. As a consequence, the water quality is closely monitored, especially regarding dangerous substances. Water analysis for dangerous substances is carried out every three months.

The dangerous substances analysed are the following: ammonium, cyanide, arsenic, cadmium, chromium, copper, nickel, lead, zinc, phenol index, total hydrocarbons, BTEX (4), PAHs (16), formaldehydes, isocyanates, PCBs, halogenated VOCs.

Mine water quality results

The treatment scheme of La Houve is located in Creutzwald, it treats mine water from the Western mine water reservoir. The treatment scheme has been operating since November 2009. In the first two years, cyanide and a few PAHs and BTEX were found in raw mine water at concentrations higher than the environmental quality standards (EQS). Since then, no dangerous substances have been quantified at concentrations higher than the EQS.

The treatment scheme of Simon treats mine water from the Central-Eastern mine water reservoir. The scheme is located in Forbach and has been operating since November 2012. Copper and zinc have been found in raw mine water at concentrations higher than the EQS, since the beginning of pumping. The other substances have not been quantified at concentrations higher than the EQS.

The treatment scheme of Vouters also treats water from the Central-Eastern mine water reservoir, like the treatment scheme of Simon. The scheme is located in Freyming-Merlebach and has been operating since July 2015. Copper and zinc have also been found in raw mine water at concentrations higher than the EQS, since the beginning of pumping. The other substances have not been quantified at concentrations higher than the EQS.

Content evolution of iron and manganese

While pumping stations enable to keep mine flooding under control, treatment schemes enable to lower the raw mine water content in iron, manganese and total suspended solids (TSS), before discharging the treated mine water in local rivers. To design water treatment schemes, predictions of mine water content are essential.

Mine water content predictions

Predictions of iron and manganese evolution in the raw mine water were made using the work of Younger [START_REF] Younger | Predicting temporal changes in total iron concentrations in groundwaters flowing from abandoned deep mines: a first approximation[END_REF]. Similar hypothesis were used for the different pumping stations and water treatment schemes.

As an example, the main hypotheses for iron evolution were:

pic concentration of 98 mg/L (± 14 mg/L) and stabilization at 5 mg/L (± 1,25 mg/L) -kinetics based on a mine water emergence with a constant flow rate -10 % of the volume of the mine reservoir to be leached out by water circulation

Predictions versus reality

Iron and manganese evolutions at the pumping station of La Houve are presented Figure 3. The predictions for iron and manganese evolution at the pumping stations turned out to be wrong. Concentrations were supposed to be high at the start of the emergence, but decreasing quite strongly in a short time. Looking at the real data, it appears that the concentrations did not go that high, and then decreased in a much slower way. Such a difference between the predictions and the real data is a major issue, as the predictions were used to design the water treatment schemes.

However, predicting iron and manganese evolution is not an easy task. Furthermore, at the time CdF did that study, the pumping stations were not yet an option, and an emergence was supposed to be set up instead. Using the work of Zeman [START_REF] Zeman | Dynamic model of long term geochemical evolution of mine water after mine closure and flooding[END_REF], predictions for iron and manganese evolution were updated by the BRGM (Figure 4).

While the initial study predicted that raw mine water would reach an iron concentration above 10 mg/L by 2013, the new prediction clearly shows that even in 2020, the iron concentration will not have reached the 10 mg/L threshold.

Mine water content predictions feedback

The main issue with the predictions made using Younger's work was the use of a constant flow rate. Indeed, the mine water pumping stations run at a variable pumping rate that keeps increasing over the years and will reach its maximum at around 2050. As shown in Figure 5, iron concentration tends to increase with the flow rate. Once the flow rate is stabilized, the iron concentration decreases.

Conclusions

Looking at the example of the French coal mine flooding, it appears that for mine flooding predictions, collecting and computing accurate data is the key element to keep in mind. The use of an isolated exploitation field as a case study to determine coefficient for residual voids, has also been decisive for the accuracy of the predictions.

Regarding the water quality of the two mine water reservoirs of the Lorraine basin, only a few toxic substances were quantified above the EQS at the start of the La Houve treatment scheme. Since then, no dangerous substances have been quantified at concentrations higher than the EQS, whether in the Western mine water reservoir, or in the Central-Eastern mine water reservoir.

At last, predicting water quality evolution in terms of iron and manganese is known to be a difficult task. In the Lorraine basin, the predictions did not take into account the evolution of the pumping rates of the pumping stations. As a consequence, the predictions overestimated the pic concentration, and underestimated the time for the concentration to decrease. Such a difference between the predictions and the real data is a major issue, as predictions of water quality are used to design water treatment schemes.

Figure 1 -

 1 Figure 1 -Mine workings and mine water reservoirs of the Lorraine-Saar basin.

Figure 2 -

 2 Figure 2 -Mine flooding of the Forbach exploitation field.

Figure 3 -

 3 Figure 3 -Iron and manganese evolution at the pumping station of La Houve.

Figure 4 -

 4 Figure 4 -Updated iron prediction at the pumping station of La Houve.

Figure 5 -

 5 Figure 5 -Correlation between variations in flow rate and iron concentration at La Houve pumping station.


Cécile Chenot
email: cecile.chenot@cea.fr

Jérôme Bobin
email: jerome.bobin@cea.fr

Blind separation of sparse sources in the presence of outliers

Keywords: Blind source separation, sparse signal modeling, robust recovery, outliers

Blind Source Separation (BSS) plays a key role to analyze multichannel data since it aims at recovering unknown underlying elementary sources from observed linear mixtures in an unsupervised way. In a large number of applications, multichannel measurements contain corrupted entries, which are highly detrimental for most BSS techniques. In this article, we introduce a new robust BSS technique coined robust Adaptive Morphological Component Analysis (rAMCA). Based on sparse signal modeling, it makes profit of an alternate reweighting minimization technique that yields a robust estimation of the sources and the mixing matrix simultaneously with the removal of the spurious outliers. Numerical experiments are provided that illustrate the robustness of this new algorithm with respect to aberrant outliers on a wide range of blind separation instances. In contrast to current robust BSS methods, the rAMCA algorithm is shown to perform very well when the number of observations is close or equal to the number of sources.

Introduction

The fast development of multichannel sensors and imagers in a wide range of scientific fields mandates the development of dedicated data analysis tools.

In this context, Blind Source Separation (BSS) plays a key role as it allows extracting relevant information in an unsupervised manner. Moreover, BSS has demonstrated its efficiency in numerous applications such as astrophysics [START_REF] Bobin | Joint Planck and WMAP CMB map reconstruction[END_REF] or hyperspectral unmixing [START_REF] Fevotte | Nonlinear Hyperspectral Unmixing With Robust Nonnegative Matrix Factorization, Image Processing[END_REF] to only name two. In the setting of BSS, the data are made of m multichannel observations {X i } i=1..m . Each one is assumed to be composed of the linear mixture of n ≤ m sources {S j } j=1..n with t > m samples. This so-called linear mixture model is generally recast in the following matrix formulation:

X = AS,
where X ∈ R m×t stands for the observations, A ∈ R m×n the mixing matrix, and S ∈ R n×t the sources. The objective of BSS is to estimate both A and S from the knowledge of the observations X only. Without any further 25 assumption, BSS is a challenging matrix factorization problem that admits an infinite number of solutions. Prior information on the sources and/or the mixing matrix is required to tackle such an ill-posed problem. In brief, most BSS methods mainly differ from the prior used to describe the sources such as statistical independence in ICA (Independent Component Analysis), 30 non-negativity in NMF (Nonnegative Matrix Factorization) or sparsity. For more details about standard BSS methods, we refer to the seminal book [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF] and references therein.

It is well-known that the presence of noise hampers the performances of most BSS methods [START_REF] Gadhok | Rotation sensitivity of independent component analysis to outliers[END_REF][START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF]. Accounting for additional additive Gaussian noise is carried out by adding an extra noise term N to the linear mixture model: X = AS+N. However, in a large number of applications, the observations are also contaminated with aberrant entries, rare and large errors, which are not correctly modeled by Gaussian noise models. More precisely, such spurious outliers include observed unexpected physical events or malfunctions of captors.

Important examples are: i) stripping noise or impulse noise in hyperspectral data [START_REF] Li | Denoising of Hyperspectral Images Employing Two-Phase Matrix Decomposition, Selected Topics in Applied Earth Observations and Remote Sensing[END_REF]), ii) cosmic ray contamination in astronomical images [START_REF] Van Dokkum | Cosmic-ray rejection by laplacian edge detection[END_REF], iii) point sources emissions in astrophysical data [START_REF] Sureau | Sparse pointsource removal for full-sky cmb experiments: application to wmap 9-year data[END_REF] to only name a few. Beyond instrumental or physical artifacts, it has been recently advocated that sparse deviations from the linear mixture model can be approximated with outliers models in hyperspectral data [START_REF] Fevotte | Nonlinear Hyperspectral Unmixing With Robust Nonnegative Matrix Factorization, Image Processing[END_REF]. Accounting for both Gaussian noise and outliers, we will further consider that the observations can be expressed as follows:

X = AS + O + N, (1)
where O ∈ R m×t stands for the outliers, and N ∈ R m×t the Gaussian noise.

Extending the BSS framework to further dealing with outliers refers to robust BSS.

Robust BSS methods in the literature

In the current literature, only few BSS methods have been developed to manage the presence of outliers. These techniques can be split into three different groups:

• Robust Independent Component Analysis : In this framework, the sources are assumed to be mutually independent. One way to measure statistical independence of the estimated sources is to compute their mutual information, defined as the Kullback-Leibler (KL) divergence between the product of their marginal distributions and their joint distribution [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]. Unfortunately, the KL-divergence is highly sensitive to the presence of outliers [START_REF] Gadhok | Rotation sensitivity of independent component analysis to outliers[END_REF]. To overcome this problem, the KL divergence is substituted in [START_REF] Mihoko | Robust blind source separation by beta divergence[END_REF] with the more robust β-divergence. While such methods provide a robust estimation of the mixing matrix, it however does not perform any sources/outliers separation.

• Two-steps methods : These methods are comprised of two successive steps: i) removal of the outliers from the observations in a first step, and ii) separation of the sources from the "outliers-cleaned"data using standard non-robust BSS techniques. The first outliers denoising step is however complex to tackle. Recently, and to the best of our knowledge, the most powerful techniques are based on the PCP algorithm (Principle Component Pursuit - [START_REF] Candès | Robust principal component analysis?[END_REF]), which have been proposed in [START_REF] Li | Denoising of Hyperspectral Images Employing Two-Phase Matrix Decomposition, Selected Topics in Applied Earth Observations and Remote Sensing[END_REF][START_REF] Zhang | Hyperspectral Image Restoration Using Low-Rank Matrix Recovery[END_REF]. However, it is essential for the PCP algorithm to work that the contribution of the sources AS has low rank. This assumption is valid whenever m n. If this assumption holds true for hyperspectral data, it is far from being the case in more general BSS problems such as in astrophysics [START_REF] Bobin | Joint Planck and WMAP CMB map reconstruction[END_REF].

• Component separation methods : These approaches aim at recovering simultaneously A, S and O. They have been essentially used in the NMF framework with optimization methods [START_REF] Fevotte | Nonlinear Hyperspectral Unmixing With Robust Nonnegative Matrix Factorization, Image Processing[END_REF][START_REF] Zhang | Robust non-negative matrix factorization[END_REF][START_REF] Shen | Robust nonnegative matrix factorization via 1 norm regularization[END_REF][START_REF] Li | Hyperspectral unmixing with robust collaborative sparse regression[END_REF][START_REF] Halimi | Fast hyperspectral unmixing in presence of nonlinearity or mismodelling effects[END_REF] or a Bayesian approach [START_REF] Altmann | Robust Linear Spectral Unmixing Using Anomaly Detection[END_REF]. These methods strongly rely on the positivity of both the sources and the mixing matrix, which is not necessarily a valid assumption in general settings. This is especially the case in imaging where the sources are more conveniently modeled in transformed domains where the non-negativity assumption does not hold.

Recently, we introduced a component separation method coined robust Generalized Morphological Component Analysis (rGMCA - [START_REF] Chenot | Robust Sparse Blind Source Separation[END_REF]), which does not require assuming that the sources and/or the mixing matrix are non-negative.

This algorithm emphasizes on the sparse modeling of the sources and outliers in the same signal representation. We showed that the rGMCA algorithm provides good separation performances in the over-determined setting (m > n) but fails at precisely solving robust BSS problems when the number of observations is close to the number of sources (determined case, m = n). This highly limits its suitability in applications where the number of available observations is of the order of the number of sources, such as in astrophysics [START_REF] Bobin | Joint Planck and WMAP CMB map reconstruction[END_REF].

Contribution:

In this article, we propose a novel algorithm, coined robust Adaptive Morphological Component Analysis (rAMCA), that generalizes the rGMCA algorithm [START_REF] Chenot | Robust Sparse Blind Source Separation[END_REF]. It first builds upon the sparse modeling of both the outliers and the sources in the same transformed domain. Unlike the rGMCA algorithm, the proposed algorithm further relies on two novel elements: i) a refined modeling of the outliers in the source domain based on an analogy with partially correlated sources (see [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF]), and ii) an improved outliers estimation procedure, described in Section 2 and 3 respectively. In Section 4, the proposed method is shown to yield a highly effective estimation of the mixing matrix. It also performs very well when the number of observations 105 is close or equal to the number of sources; a challenging setting for which currently available robust BSS methods fail. Besides, we describe how the parameters of rAMCA can be automatically tuned.

Notations

Uppercase boldface letters denote matrices. The Moore-Penrose pseudoinverse of the matrix M is designated by M † . The jth column of M is denoted M j , the ith row M i , and the i, jth entry M i,j . The norm M 2 denotes the Frobenius norm of M, and more generally M p designates the p-norm of the matrix M seen as a long vector. The soft-thresholding operator is denoted S λ (M), where

[S λ (M)] i,j =    M i,j -sign(M i,j ) * λ i if |M i,j | > λ i 0 otherwise
Last, the operator MAD designates the median absolute deviation and Pr 110 stands for probability.

Sparse BSS in the presence of outliers

Impact of outliers on sparse BSS methods

Based on sparse modeling [START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF], sparse BSS assumes that the sources {S} i=1..n have sparsely distributed entries in some signal representation Φ:

S i = α i Φ, ∀i = 1..n,
where α i is composed of a small number of non-zero elements (i.e. exactly sparse model) or few significant large-amplitude entries (i.e. approximately sparse model). Most natural signals verify such a sparsity property in a wellsuited signal representation Φ such as the wavelets or the curvelets to only name two [START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF]. Sparsity has been shown to largely improve the performances of BSS methods [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF][START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF][START_REF] Zibulevsky | Blind Source Separation by Sparse Decomposition in a Signal Dictionary[END_REF] since it allows for an enhanced discrimination between the sources to be estimated, which are assumed to verify some morphological diversity principle (MDP) [START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF]. According to the MDP, the largeamplitude coefficients of the sources in Φ, which encode their most salient morphological features, are distinct. For instance, sparse and statistically independent sources verify the MDP with high probability. This is illustrated in fig. 1b: the source's samples in Φ (blue dots) with the largest amplitudes are mainly clustered along the canonical axes. Based on this principle, the Generalized Morphological Component Analysis (GMCA) algorithm [START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF] tackles sparse BSS problems by seeking the sources that are jointly the sparsest possible. This is performed by minimizing:

minimize S,A 1 2 X -AS 2 2 + n i=1 λ i S i Φ T 1 ,
where the first term is the data-fidelity term and the second term promotes the sparsity of the sources in the transformed domain Φ.

In the following, for the sake of clarity, we will assume that the sources and the outliers are sparse in the direct domain Φ = I. Similarly to [START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF] and [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF],

all the results will be exact for any orthogonal transform Φ and a good approximation for redundant sparse representations such as tight frames, which have diagonally dominant Gram matrices (e.g. undecimated wavelet transforms, curvelets [START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF], etc.).

Unlike pure sources' samples, the outliers are assumed to be distributed in general position: they do not cluster in any specific direction as illustrated in fig. 1a (red dots). Since the outliers are sparse and distributed in general position, we will further assume that few columns of O are entirely active such as in [START_REF] Fevotte | Nonlinear Hyperspectral Unmixing With Robust Nonnegative Matrix Factorization, Image Processing[END_REF].

According to the MDP, each of the sources is precisely described by its most prominent coefficients, which makes them the most informative samples to tackle the separation task. However, in the presence of outliers, the largeamplitude entries of the projected sources Ŝ = A † X are more likely related to samples that are corrupted by outliers as shown in fig. 1b. Consequently, when seeking for the jointly sparsest sources by unmixing the large entries of X, sparse BSS methods are misled by the outliers. This is illustrated in the fig. 1c which shows the scatter plot of the estimated projected sources S = à † X where à is the mixing matrix returned by GMCA: the corresponding sources are jointly sparser than the ones in fig. 1b. However such a solution yields a poorly estimated mixing matrix Ã. It highlights that applying sparse BSS method will very likely yield an erroneous solution. By nature, sparse BSS methods are highly sensitive to the presence of outliers.

An analogy with partially correlated sources

In [START_REF] Chenot | Robust Sparse Blind Source Separation[END_REF], it clearly appears that solving robust BSS problems requires discriminating between the outliers and the sources. Since both the outliers and the sources are assumed to be sparsely represented in the same domain If the mixing matrix A were perfectly known, the sources would be approximated by projecting the corrupted data onto the span of A: Ŝ = A † X. The estimated sources are the linear combination of the clean sources and the projected outliers: Ŝ = S + Ô. Due to the projected outliers contribution which is broadly distributed, some of the largest entries of Ŝ are active simultaneously in several sources (c.f. the red contribution in fig. 1b). These shared active samples are reminiscent of the partial correlations of the sources we discussed in [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF]. Indeed, the samples of partially correlated sources can be similarly divided into two groups: the discriminant samples respecting the MDP (the jointly sparse contribution in blue in fig. 1b) and the samples corresponding to the partial correlations which active simultaneously in several sources (the broadly distributed contribution in red in fig. 1b). Unlike the rGMCA algorithm we introduced in [START_REF] Chenot | Robust Sparse Blind Source Separation[END_REF], we propose to exploit the analogy between the impact of the projected outliers and sparse and partially correlated sources, which yields a novel robust BSS algorithm that is described in the next section.

Robust AMCA Algorithm

Following the analogy between the impact of outliers and partial correlations, the rAMCA algorithm will be built upon the AMCA algorithm (Adaptive Morphological Component Analysis), which has been designed to deal with partially correlated sources [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF]. In this specific context, we underlined in [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF] that the ability to identify correlated entries is critical to perform the separation. For that purpose, the AMCA algorithm builds upon an adaptive weighting scheme that assigns to each column of the observation coefficients X a weight, whose goal is to penalize correlated entries in the separation process. Details about the weighting procedure will be given below. According to [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF], in the setting of partially correlated sources, the AMCA algorithm performs by minimizing the following problem: minimize

A,S 1 2 (X -AS) W 2 2 + n i=1 λ i S i 1 , (2)
where W ∈ R t×t is the weight matrix.

In the spirit of [START_REF] Chenot | Robust Sparse Blind Source Separation[END_REF], we propose to estimate jointly A, S and O by exploiting the sparsity of the sources and the outliers. Unlike the rGMCA algorithm, we propose to further employ a weighting scheme similar to AMCA. This can be done by substituting the problem in eq. 2 with the following one: minimize

S,A,O 1 2 (X -AS -O) W 2 2 + n i=1 λ i S i 1 + β O 2,1 . (3)
The 2,1 norm, defined such as

O 2,1 = m j=1
O j 2 , favors solutions O with few entirely active columns. This regularization term is well suited to capture outliers that are distributed in general position in the data domain.

This problem is non-convex but can be tackled using a minimization procedure such as the Block Coordinate Relaxation [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF] (BCR). This minimization technique amounts to sequentially minimizing subsets of variables. A natural choice would consist in estimating alternatingly the three variables of interest A, S, and O. However, we found that this choice performs poorly in practice since errors are more likely to propagate from one variable to the other during the sequence of minimization steps. We rather opted for sequential minimization of two blocks of variables: i) the couple (A, S) and ii) the outliers matrix O. The major advantage of this choice is that it provides a much more robust estimation of (A, S) since both parameters are updated jointly and not independently from the residual X -O.

Procedure 1 rAMCA Algorithm

1: procedure rAMCA(X, n) 2:
Initialize Ã(0) (randomly or with a PCA), S(0) = 0 and Õ(0) = 0. Update S(i,k) from ( 5)

7:
Update W from ( 7)

8:
Update Ã(i,k) from ( 6)

9: Set S(k) ← S(i-1,k) and Ã(k) ← Ã(i-1,k) 10:
Update Õ(k) from ( 9) Estimation of O return S(k-1) , Ã(k-1) , Õ(k-1) . Applying the BCR technique to estimate the mixing matrix and the sources amounts to minimizing the problem in Eq. 3 assuming O is fixed: minimize

Estimating the sources and the mixing matrix

S,A 1 2 (X -AS -O) W 2 2 + n i=1 λ i S i 1 . (4)
The problem shares similarities with the problem solved by the AMCA algorithm (see [START_REF] Fevotte | Nonlinear Hyperspectral Unmixing With Robust Nonnegative Matrix Factorization, Image Processing[END_REF]) with the exception that it applies to the residual X -O rather than the raw observations X. Following the AMCA algorithm, the 205 problem in ( 4) is tackled by minimizing alternately the cost function with respect to A and S with the two following steps:

• Updating S assuming A is fixed : Minimizing (4) with respect to S consists in solving the following convex problem: minimize

S 1 2 (X -AS -O) W 2 2 + n i=1 λ i S i 1 .
Unless A is orthogonal, the previous problem does not admit a closed form solution. In the spirit of alternated least-square minimization techniques, we proposed in [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF] to rather approximate this step with a projected least-square, which highly limits the computational cost of the update:

S i = S λ i A † (X -O) i . (5)
• Updating A assuming S is fixed : Minimizing (4) with respect to A amounts to solving the following convex problem: minimize

A 1 2 (X -AS -O) W 2 2 .
which admits a closed form solution:

A = (X -O) W (SW) † . (6)
In practice, to avoid the balance indeterminacy between A and S, we assume that the columns of A are normalized for the 2 norm. Similarly to what is done with AMCA, this additional constraint is handled by normalizing the columns of A after the projected least-squares.

Similarly to the AMCA algorithm [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF], the weights play a central role. In the setting of robust BSS, they help providing robustness to the remaining outliers contribution in the estimated residual X -O. Following the analogy with partial correlations, the weights aim at penalizing entries of the estimated sources which are in general position rather than clustered along a canonical axis. The former are more likely related to residuals of outliers while the latter are characteristics of the sources. Following [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF], samples in general position can be traced by measuring the sparsity level of the columns 225 of the estimated sources using a q norm. Therefore, the diagonal elements of the weight matrix W are defined as follows:

W t,t = 1 S t q + , (7)
where S denotes the normalized sources

S i = S i S i 2
and where is a scalar typically small used to avoid numerical issues. The parameter q is chosen in the range [0, 1]. For more details, we refer the reader to [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF].

O 1 2 (X -AS -O) W 2 2 + β O 2,1 .
Given that only the diagonal terms of W are non-zero, this problem is separable. It amounts to solve for each sample k ∈ {1..t}:

minimize O k 1 2 (X -AS) k -O k W k k 2 2 + β O k 2 .
This problem is equivalent to:

minimize O k (W k k ) 2 2 (X -AS) k -O k 2 2 + β O k 2 .
Then, by setting β = β (W k k ) 2 , we end up with:

minimize O k 1 2 (X -AS) k -O k 2 2 + β O k 2 .
This problem has a closed form solution which has been derived in [START_REF] Kowalski | Sparse regression using mixed norms[END_REF]:

O k = (X -AS) k × 1 - βk (X -AS) k 2 + . ( 8
)
Detecting the outliers. Most sparsity-based thresholding procedures can be interpreted as detection procedures: detecting sparse samples out of dense noise. In that case, it is customary to fix the value of the threshold based on the noise statistics [START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF]; we will see later that this is exactly how the thresh- In that case, only the Gaussian noise contributes to the residual outside the support of O. Therefore, the samples (X -AS) k

olds {λ i } i=1,••• ,
2 k: O k 2 =0
follow a χ law with m degrees of freedom. The value of β can then be chosen based on the expected value of the χ law:

σ × √ 2 × Γ( m+1 2 ) Γ( m 2 )
, where σ corresponds to the standard deviation of N.

Nevertheless, relying on the noise statistics only provides a detection procedure that is not reliable in the determined case. Indeed, even if A is correctly recovered, the outliers are very likely to leak into the estimated sources S since they also lie in the span of A: S = S + A † O, such that

A S = AS + O.
An accurate detection of the outliers based on the residual X -A S is then not possible. To overcome this issue, we propose to rather build the detection procedure on a quantity that allows discriminating between the outliers and the sources, especially in the determined case.

We emphasized in Section 2.2 that in the source domain the entries of S are jointly sparse, i.e. clustered along the canonical axes, whereas the projected outliers behave as correlated non-sparse entries. In this context, the δ-density, which has been introduced in [START_REF] Studer | Recovery of Signals with Low Density[END_REF], provides a convenient measure of sample sparsity that allows discriminating between sparse and non-sparse columns of S. The δ-density of any jth non-zero sample of the estimated sources is defined as

δ( Sj ) = Sj 1 Sj ∞
. This quantity takes its values between 1 (for one active entry) and n (for a column whose entries have the same amplitudes). More interestingly, it is independent of the amplitude of the columns and well suited for sparse and approximately sparse signals. In this framework, detecting the support of O can be performed by identifying the columns of the estimated S whose δ-density is larger than a certain threshold α that needs to be determined. This is somehow reminiscent of the outliers detection discussed in [START_REF] Wright | Robust face recognition via sparse representation[END_REF].

In the general setting, determining an optimal numerical value for α is challenging without an accurate statistical modeling of the sources and the outliers. In the following, we propose to use the following statistical modeling:

• the sources are drawn from a Generalized Gaussian law with parameter ρ denoted by G(ρ).

• the amplitude of the outliers in the sources domain follows a Gaussian law N , well suited to model samples that are distributed in general position.

Let us notice that the variances of these statistical models do not matter since the δ-density is independent of the amplitude. From this statistical model, the threshold α is derived from a classical hypothesis testing procedure such that, for any random variable X of size n:

Pr (δ(X) < α| X ∼ G(ρ)) = Pr (δ(X) > α| X ∼ N ) .
where Pr (δ(X) < α| X ∼ G(ρ)) stands for the probability for the δ-density to be smaller than α assuming that every entry of X is distributed according to a Generalized Gaussian law with parameter ρ. According to [START_REF] Sureau | Sparse pointsource removal for full-sky cmb experiments: application to wmap 9-year data[END_REF], the amplitude of the detected outliers is derived from the estimated residual X-ÃS . Previously, we underline that X-ÃS is very likely to contain some errors. A more conservative but more effective choice consists in deriving the amplitude of the detected outliers from the data X.

As a summary, the outliers O are estimated as follows:

Õk =        0 if δ( Sk ) < α X k × 1 - β X k 2 otherwise, (9)
where

β = MAD(X -AS -O) × √ 2 × Γ( m+1 2 ) Γ( m 2 )
and MAD(X -AS -O) corresponds to a good estimate of the standard deviation of N if it is not known.

Despite the simplicity of the statistical model used to derive a value for α and consequently β, the proposed scheme has been proved to be robust in the various numerical experiments of Section 4. Furthermore, at each iteration of the rAMCA algorithm 1, the couple (A, S) is fully re-estimated, which also makes the algorithm less sensitive to mis-estimations of the outliers O.

Choice of the parameters

Strategy for λ: The major parameters of the sparse source separation problems are the thresholds λi i=1..n . Similar to [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF], we use the decreasing thresholding strategy proposed in [START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF] which has two interesting properties: i) it prevents the incorporation of noise in the source estimates and ii) it makes AMCA less prone to be being trapped in local minima.

Robustness to Gaussian noise: The soft-thresholding operator S λ (.) rejects the entries having an amplitude smaller than λ. The final threshold is thus chosen based on the level of noise which contaminates the projected sources A † (X -O) = S + A † N, so as to remove the additional noise. Indeed, Gaussian noise removal (A † N) from sparse signals (S) based on thresholding can be interpreted as a standard hypothesis testing [START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF]. The value of λ i is typically set to 3σ i , where σ i stands for the standard deviation of the noise contaminating the ith projected source, i.e. (A † N) i . If these values are not known, they can be estimated empirically using the Median Absolute Deviation (MAD) since MAD A † (X -O) i ≈ MAD (A † N) i for sparse sources -see [START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF][START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF].

Robustness to local minima: Following [START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF], the use of a decreasing thresholding strategy remarkably improves the separation performances since it provides more robustness to the spurious local minima. During the unmixing process, the thresholds are chosen automatically so that the number of non-zero entries of the sources is increased by a fixed amount at every iteration. More precisely, given the total number of iterations I, the jth projected source S(i,k)

j = Ã(i-1,k) † (X -Õ(k) )
j is thresholded at the ith iteration by:

λ j = pct | S(i,k) j | | S(i,k) j |>3σ j , 100 × I -i I ,
where pct(x, v) denotes the vth percentile of the entries of x.

Number of inner loops I:

The number of iterations is set to I = 1000, which turned to be a good compromise in the numerical experiments.

Strategy for β: In the spirit of the decreasing value strategy used for λ in AMCA, the number of eligible active samples of the estimated outliers is increased during the algorithm. More precisely, at the kth iteration, we select the outliers among the 5k% largest entries of the residue in order to limit the number of false estimations. We underline that these parameters are also automatically determined: α depends only on the number of sources and β on the number of observations.

Number of outer loops K:

Last, the number of outer loops is maximally set to 100. In practice, the algorithm is stopped when à and Õ are jointly stabilized fig. 3. More precisely, rAMCA stops at the kth outer loop if: max j=1..n Ã(k-1)j , Ã(k)j < 5°, and supp Õ(k-1) = supp Õ(k) , where supp(x) denotes the support of the vector x.

Stability of rAMCA

Since the problem (3) is not convex, we can only expect to converge to a local minima. Besides, given that the proposed strategy uses varying parameters, the convergence to a critical point, strictly speaking, cannot be proved.

However, the stability of the two variables of interest, the support of the corrupted samples and the mixing matrix, is heuristically well motivated.

We propose to minimize the function using the Block Coordinate Descent (BCD) method [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]. It has been shown in [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF] that minimizing (3) alternately for each variable with fixed parameters converges to a stationary point. However, in practice, minimizing (3) with the cyclic rule and with fixed parameters performs poorly: this minimizing scheme is very likely to be prone to being trapped in local minima. That is why, we minimize 3 using a sequential minimization alternating between the blocks (A, S) and O, as well as the decreasing parameters strategy.

Once the detrimental outliers (or the data estimated as being detrimental)

have been removed from the observations, the AMCA algorithm, whose stability has been discussed in [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF], returns a similar A from one iteration to another (since the input X -O is constant from one outer iteration to another one), fig. 3a,3b.

For illustrative purpose, we display the maximal angle made between the columns of Ã(k) and Ã(k+1) (see Section 4 for the metrics) as well as the percentage of estimated corrupting columns for n = m = 10 sources generated according to Section 4.1 and 30% of corrupting columns. After few outer loops, the number of estimated columns fig. 3a andà almost not vary fig. 3b (a variation with the maximal order of magnitude of 10 -3 is observed for A due to the projected least squares).

Numerical Experiments

Experimental protocol

In this section, rAMCA is compared with various robust BSS algorithms: • GMCA [START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF]: this standard sparse BSS method is used to illustrate the sensitivity of the non-robust BSS algorithms to the presence of outliers.

• AMCA [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF] whose performances show the benefits of the weighting scheme (difference between AMCA and GMCA) and of the explicit 350 estimation of O (difference between AMCA and rAMCA).

• rGMCA [START_REF] Chenot | Robust Sparse Blind Source Separation[END_REF]: the discrepancy between its performances and the ones of rAMCA illustrates the key roles of the novel penalization and outliers detection procedure, which are, unlike rGMCA, based on the refined modeling of the outliers in the source domain.

355

• the robust minimization of the β-divergence [START_REF] Mihoko | Robust blind source separation by beta divergence[END_REF], (implementation similar to [START_REF] Gadhok | An Implementation of β-Divergence for Blind Source Separation[END_REF]), which assumes that m = n and only estimates the mixing matrix.

• the robust combination PCP+GMCA: the outliers are first estimated with PCP [START_REF] Candès | Robust principal component analysis?[END_REF] which assumes that m n, and then the sources and mixing matrix are estimated with GMCA.

The parameters of PCP+GMCA and of the minimization of the β-divergence are manually tuned. In the first part of this section, their performances are evaluated on various scenarios with synthetic data, which allows performing Monte-Carlo simulations.

Complexity of the different methods. The AMCA and GMCA based methods have a complexity of O (mnt), while PCP has a complexity given by O (m 2 t).

Their complexities are thus similar except if m n, what is required by PCP (low-rank assumption). The minimization of the β-divergence depends on the algorithm used to perform the minimization.

Nonetheless, we point out that in practice, if the dimensions are moderate, running rAMCA may require more computational time than running once the combination PCP+GMCA. On the other hand, PCP is a parametric method whose parameter tuning requires several trials, what is then even more time consuming.

Performance criteria. We emphasized in [START_REF] Chenot | Robust Sparse Blind Source Separation[END_REF] that the algorithms listed above do not all yield a precise estimation of the sources but rather provide a robust estimation of A. Therefore, we will focus on assessing the performances of these algorithms with respect to the mixing matrix. More precisely, we propose to evaluate the accuracy of the different algorithms as well as their reliability, which is particularly relevant since BSS problems are non-convex.

The quantity ∆

A = P Ã † A-I 1 n 2
is used as a global indicator of the mixing ma-trix estimation accuracy [START_REF] Bobin | Sparsity and morphological diversity in blind source separation, Image Processing[END_REF], where the matrix P corrects for the permutation indeterminacy. Additionally, for every simulation and for each algorithm, we record the number of runs for which A has been correctly recovered (normalized to 1). The mixing matrix is said to be correctly recovered if, for every column of A, the angle between the estimated and true ith column is smaller than 5 • : arccos( Ãi , A i ) < 5°. This quantity provides a good indicator of the reliability of the algorithms.

Data Setting. The comparisons are first carried out on synthetic data in order to illustrate the impact of parameters such as the percentage of corrupted data or the number of observations with Monte Carlo simulations (48 simulations). The data are generated as follows:

• A total of 8 sources (unless otherwise stated) are drawn from a Bernoulli-Gaussian law whose activation rate is fixed to 5%, and the standard deviation of their amplitude σ S to 100. The number of samples t is fixed to 4096.

• The mixing matrix is drawn according to a normal law with zero mean.

The columns of A are normalized to unit 2 norm.

• The outliers are generated so as to corrupt at random a low number of columns of X. The activation of these columns is drawn according to a Bernoulli process with probability ρ, which fixes the average number of corrupted columns to ρt. The amplitude of the outliers is drawn at random from a Gaussian distribution with zero mean and standard deviation σ O .

• The noise is generated according to a Gaussian distribution with zero mean. Its standard deviation is set to 0.1.

Influence of the number of observations

We emphasized in [START_REF] Chenot | Robust Sparse Blind Source Separation[END_REF] that the separation of the sources contribution and the outliers is more challenging if m is close to n. The ratio m n is therefore a crucial parameter in BSS, especially in the presence of outliers. In this paragraph, the data are composed of m observations. The amplitude of the outliers is fixed to σ O = 100 for n = m and then the amplitude ratio between the outliers and the sources contribution is kept constant. The percentage of outliers is fixed to 10% with ρ = 0.1.

As shown in fig. 4a, rAMCA tends to be less influenced by the number of observations. The results of all the methods (except the β-divergence minimization algorithm) are better if m is very large: the condition number of A is smaller and the outliers can be better distinguished from the sources contribution since the energy of the outliers lying in the subspace generated by A is lower when m is large. In this regime, the low-rankness of the term AS becomes a valid assumption, which makes PCP more efficient [START_REF] Candès | Robust principal component analysis?[END_REF].

The results are not strictly improved with an increasing number of measurements for the β-divergence algorithm. Since the β-divergence minimization algorithm has been designed for the determined case only, its application to the over-determined case requires a first dimension reduction step. This pre-processing step, which is performed by PCA, is also impacted by the presence of outliers and hampers the performances of this algorithm.

In order to further illustrate the impact of the ratio m n , the errors As shown in fig. 5, all the algorithms are less reliable if only few samples are available since the clustering aspect of the sources contributions is not significant (an unmixing, even without outliers, is challenging if only few samples are available). Besides, all the strategies become more and more precise as the number of samples t increased.

Increasing the number of samples has several favorable effects on the unmixing: the number of samples available to unmix the sources becomes sufficient regardless of the presence of outliers, and the clustered aspect of AS has a greater importance since there are more and more clustered samples in the term AS but the outliers are still in general position (generating randomly several outliers in a same direction is quite unlikely).

The results would have been different if the proportion of data samples in a given direction had been set constant from one value of t to another. For instance, if one resizes an multi/hyperspectral data cube, these proportions are kept constant, and for the largest image size, few but several outliers are in a same direction. There are some applications (most of the observations of physical processes) were the outliers are not strictly speaking in general position (such as in 4.6), but whose contributions are less structured/clustered than the one of the sources: the weighting scheme penalizes the less clustered solutions, and so, still returns A. That is why AMCA and rAMCA requires less samples than the others methods to perform accurately fig. 5.

In the following, the impact of two other parameters will be investigated: the percentage of corrupted data and their amplitudes. We will focus on the determined case which is more challenging. Since the low-rankness assumption makes no sense in the determined case, the algorithm PCP+GMCA will not be evaluated.

Influence of the amplitude of the outliers

In the following experiments, we consider that 10% of the data samples are corrupted with outliers. Fig. 6a shows the behavior the algorithms when the amplitude of the outliers σ O varies.

The figure 6a shows that the standard GMCA rapidly fails to correctly recover the mixing matrix when the amplitude of the outliers increases. In these experiments, the algorithms AMCA and β-divergence minimization algorithms provide very similar results. Interestingly, rAMCA tends to be the least impacted by the amplitude of the outliers, especially when their amplitude is of the order of the source's level or very large. When the amplitude of the outliers and the sources are close, the weighting schemes of AMCA and rGMCA are less effective at penalizing the outliers. Unlike AMCA, the rAMCA algorithm progressively removes a certain level of the outliers' component, which further enhances the separation performances.

Influence of the percentage of corrupted data

In this section, the amplitudes of the outliers σ O is fixed to 100. The figure 6b shows the behavior of the BSS algorithms when the percentage of corrupted columns ρ varies according to the values of the x-axis.

As illustrated in fig. 6b, the β-divergence algorithm is able to recover correctly the mixing matrix when the number of corrupted columns of X is low (i.e. typically below 10%). The rGMCA algorithm is rapidly impacted by an do not provide satisfactory separation results. This experience is particularly challenging for these methods since: the low-rank assumption is not valid, the "sum-to-one" constraint necessary to the separation between AS and O on S for rNMF has been removed, and the outliers are less and less sparse as the width of the kernel increases. As well, the rGMCA provides good separation results when the width is low but it rapidly yields incorrect results when the width of the kernel increases. Indeed, let us recall that the outliers are also approximately sparse, which makes these separation scenarios close to the cases we investigated previously where the number of outliers is very large. This is typically the kind of settings where these methods tend to fail. The rAMCA and AMCA provide the most accurate estimates of the mixing. The discrepancy with respect to the other algorithms is particularly large when the kernel has a large width. In this regime, the level of correlation between the sources increases, a phenomenon to which the AMCA algorithm is robust [START_REF] Bobin | Sparsity and Adaptivity for the Blind Separation of Partially Correlated Sources, Signal Processing[END_REF]. Last, one of the sources estimated by rAMCA and PCP+GMCA are displayed in fig. 7b. Contrary to PCP+GMCA, the source is correctly recovered by rAMCA outside the support of O because A is correctly estimated by rAMCA. However, the leakages from the outliers into the sources estimated by rAMCA are still important: they come from the coarse scale of the wavelet coefficients, which is not sparse and for which we cannot differentiate the two contributions. Taking into account the non-negativity of the signals would limit these leakages, but necessitates the use of proximal algorithms [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] if combined with sparsity in a transformed domain [START_REF] Rapin | NMF with Sparse Regularizations in Transformed Domains[END_REF]. Nonetheless, the weighting scheme of rAMCA and AMCA is sufficient to obtain a robust estimation of A.

Software

Following the philosophy of reproducible research [START_REF] Buckheit | Wavelets and Statistics[END_REF], a python implementation of the algorithms introduced in this article will be available at http://www.cosmostat.org/GMCALab.

Conclusion

In this article, we introduce a new algorithm for tackling BSS problems in the presence of outliers, which is a key problem in a large number of ap-

(a) X 1

 1 versus X 2 . (b) Ŝ1 versus Ŝ2 .(c) S1 versus S2 .

Figure 1 :

 1 Figure 1: Three sparse sources are mixed into 4 noisy observations. Fig.(a): scatter plot of two noisy observations, (b): scatter plot of two of the projected data given by Ŝ = A † X,(c): scatter plot of the estimated sources given by S = à † X, where à has been estimated with GMCA and is far from the initial A. The initial source contribution is represented in blue, and the one of the outliers with the red stars.

4 :

 4 Set S(0,k) ← S(k-1) and Ã(0,k) ← Ã(k-1)5:while i < I do Joint estimation of A and S 6:

  200

230 3 . 2 .

 32 Estimating the outliersIn the rAMCA algorithm, the estimation of O given A and S is carried out by solving the problem in (3): minimize

  n are fixed. Similarly, and according to[START_REF] Sureau | Sparse pointsource removal for full-sky cmb experiments: application to wmap 9-year data[END_REF], the support (i.e. the set of active columns) of O is defined by the set of columns whose 2 norm exceeds the threshold β. Ideally, the columns having a 2 norm smaller than β should correspond to the remaining Gaussian noise. Consequently, the values of β should also be fixed based on the Gaussian noise statistics.

Figure 2 1 Figure 2 :

 212 Figure 2: Numerical approximations of the cumulative distribution functions of δ(.) for different values of ρ and n = 10. In blue: fig.(a): P r(δ(X) < γ|X ∼ G(1)), fig.(b) P r(δ(X) < γ|X ∼ G(0.5)), fig.(c) P r(δ(X) < γ|X ∼ G(0.3)). In red: P r(δ(X) > γ|X ∼ N ).

  Variation of A in degree across the outer loops.

Figure 3 :

 3 Figure 3: Convergence of rAMCA.

  Errors made on S and Õ versus m n .

Figure 4 :

 4 Figure 4: Influence of the number of observations on the estimations of A, S and O.

2 are

 2 displayed for a single example. Since the minimization of the β-divergence does not explicitly return O and S, we (re)-estimate O and S by minimizing (3) for fixed A, the mixing matrices estimated by the different algorithms. A good separation of S and O is possible if m n because the outliers are less likely to lie in the span of A; this is clearly shown in fig.4b. Despite an accurate recovery of A for rAMCA when m is small, the error made on the estimated outliers and sources is large fig.4b: the separation is not possible without any additional assumption on the sourcesand the outliers. Moreover, these errors decrease when the ratio m n increases whereas the error made on A remains more stable: the separation benefits from enhanced estimation of A as well as from a lower contribution of the outliers in the range of A.4.3. Influence of the number of samples tIn the following experiment, we investigate the influence of the number of samples. In order to observe the impact of this data dimension on the combination PCP+GMCA, we consider that 6 sources are mixed into 30 observations (the low-rank assumption is valid), which are corrupted by ρ O = 10% of active outliers with σ O = 50. We set σ N to 0.1. The number of samples t varies according to the x-axis of fig.5.

Figure 5 :

 5 Figure 5: Performance results of the methods versus the number of samples t.

10 0 8

 108 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 Amplitude Influence of ρ.

Figure 6 :

 6 Figure 6: Influence of the amplitude and the activation rate of the outliers.

  Influence of the kernel width.

Figure 7 :

 7 Figure 7: Top: illustration of one observation X i , without (left) and with outliers (right, corrupted entries are represented with the red dashed line). Bottom: estimated sources with rAMCA and PCP+GMCA for a width of the kernel of 6 (left) and right, performances of the different algorithms versus the width of the kernel (right).

  plications. The proposed rAMCA algorithm performs by estimating jointly the mixing matrix, the sources and the outliers. Inspired by the AMCA algorithm, it first provides a robust estimation of the sources and the mixing matrix. Additionally, it exploits the difference of structures of the outliers and the sources to provide a robust detection and estimation of the outliers based on their sparsity level in the source domain. Numerical experiments have been carried out on Monte-Carlo simulations with various experimental scenarios, which show that rAMCA yields a robust and reliable estimation of the mixing matrix. It provides the state-of-the-art separation results especially in the highly challenging determined case. Future work will exploit the difference of morphology between the sources and the outliers that is manifested in various imaging problems to further perform an accurately separation the two contributions.

Acknowledgment

This work is supported by the European Community through the grants PHySIS (contract no. 640174), DEDALE (contract no. 665044) and LENA

increasing number of corrupted data. On the other hand, the AMCA-based algorithms are less influenced by the percentage of outliers. The rAMCA algorithm provides a significantly better estimate of the mixing matrix when the number of outliers is larger than 10%.

Application to NMR spectra unmixing

In this section, we propose to compare the different algorithms in a more realistic setting: the separation of Nuclear Magnetic Resonance (NMR) spectra. In the context of spectroscopy, BSS allows to identify the different molecules of the observed mixture [START_REF] Rapin | NMF with Sparse Regularizations in Transformed Domains[END_REF]. The presence of instrumental artifacts is very frequent and makes difficult the interpretation of the data. Such artifacts can be approximated by outliers contaminating entire columns of the data matrix [START_REF] Rapin | Application of Non-negative Matrix Factorization to LC/MS data[END_REF], which is the case we investigate in the present article.

Following [START_REF] Rapin | NMF with Sparse Regularizations in Transformed Domains[END_REF], the sources are composed of 6 theoretical NMR spectra of the cholesterol, folic acid, adenosine, oleic acid, menthone and saccharose extracted from the SDBS database 1 with t = 2048 samples. These spectra are further convolved with a Laplacian kernel of varying width at half maximum (implementation from pyGMCA 2 ), which models the resolution of the instrument, fig. 7a. The set of corrupted data samples is fixed to 10 blocks of 20 consecutive columns. Their amplitudes are drawn according to a Chi-law with 1 degree of freedom, and they are further convolved with the same kernel than the sources. The amplitude of the outliers is set so that the energy of each block of outliers corresponds to the average contribution of a source in the observations Given that all the variables are non-negative, we will also compare AMCA and rAMCA with rNMF [START_REF] Fevotte | Nonlinear Hyperspectral Unmixing With Robust Nonnegative Matrix Factorization, Image Processing[END_REF], whose code is online. This method exploits the low-rankness of AS, the non-negativity and the "sum-to-one"constraint (that is, the amplitudes of each sample of S sum to one) to differentiate between the low-rank subspace and the outliers. The "sum-to-one"constraint, which is not a valid assumption in this setting, is replaced by the constraint on the columns of A, which are assumed to be normalized. We use the following inputs for rNMF: the ground truth A, the projected sources A † X + and the non-negative part of the corresponding residue.

The resulting sources admit a sparser distribution in the wavelet domain.

Subsequently, the data are transformed with the undecimated wavelet transform [START_REF] Starck | The undecimated wavelet decomposition and its reconstruction[END_REF] prior to applying the BSS algorithms, except for rNMF. Let us notice that a same wavelet transform is used for the outliers and the sources because they have a similar morphology in the present setting. In the previous experiments, we evaluated the separation performances of the algorithms in the case of exactly sparse signals. The NMR sources we consider in this section rather exhibit an approximately sparse distribution in the wavelet domain. We propose to evaluate the behavior of the robust BSS algorithms when both the sources and the outliers follow an approximate rather than exact sparse model. A simple way to evaluate the behavior of the algorithms with respect to the sparse model is to evaluate their performances when the