J Bobin

J Rapin

A Larue

J.-L Starck

Sparsity and adaptivity for the blind separation of partially correlated sources

Keywords: MLR-SSEP Source separation, sparse representations, morphological component analysis, partially correlated sources

Blind source separation (BSS) is a very popular technique to analyze multichannel data. In this context, the data are modeled as the linear combination of sources to be retrieved. For that purpose, standard BSS methods all rely on some discrimination principle, whether it is statistical independence or morphological diversity, to distinguish between the sources. However, dealing with real-world data reveals that such assumptions are rarely valid in practice: the signals of interest are more likely partially correlated, which generally hampers the performances of standard BSS methods. In this article, we introduce a novel sparsity-enforcing BSS method coined Adaptive Morphological Component Analysis (AMCA), which is designed to retrieve sparse and partially correlated sources. More precisely, it makes profit of an adaptive re-weighting scheme to favor/penalize samples based on their level of correlation. Extensive numerical experiments have been carried out which show that the proposed method is robust to the partial correlation of sources while standard BSS techniques fail. The AMCA algorithm is evaluated in the field of astrophysics for the separation of physical components from microwave data.

I. INTRODUCTION

With the rapid development of multi-wavelength or multi-sensor instruments, it is crucial to design analysis tools to efficiently extract meaningful information. This is especially the case in the field of astrophysics where the development of dedicated component separation methods has been recently of prime importance for the processing of microwave surveys such as NASA-WMAP or ESA-Planck [START_REF] Ade | Planck 2013 results. XII. Component separation[END_REF]- [START_REF] Bobin | Sparse component separation for accurate cosmic microwave background estimation[END_REF]. For that purpose, Blind source separation (BSS) is a particularly well-suited framework to analyze multichannel data. In this context and following the socalled linear mixture model, each observation or channel {x i } i=1,••• ,m is modeled as the linear combination of n ≤ m sources {s j } j=1,••• ,n :

∀i = 1, • • • , m; x i = n j=1 a ij s j + z i ,
where a ij quantifies the contribution of the j-th source in the observation x i . Each datum {x i } i=1,••• ,m and source {s j } j=1,••• ,n are assumed to have T entries. The term z i stands for additive noise or model imperfections. The well-known linear mixture model is more conveniently recast in the following matrix form:

X = AS + Z, (1)
where X is a m × T matrix the rows of which are composed by the observations, A is the so-called mixing matrix and S is the n × T matrix which contains the sources. The goal of blind source separation is to blindly estimate both the mixing matrix A and the sources S from the knowledge of X only. If one neglects the noise Z, there exists an infinite number of couples (A, S) which satisfy X = AS. This is an ill-posed problem which requires using prior information about either the mixing matrix or the sources to be recovered. For that purpose, classical approaches rely on information which promote some discriminant information or diversity between the sources. Since the beginning of the 90s, Independent Component Analysis or ICA (see and [START_REF] Comon | Handbook of blind source separation[END_REF]- [START_REF] Choi | Blind source separation and independent component analysis: a review[END_REF] references therein) has taken the lion share. In the framework of ICA, the sources are assumed to be statistically independent. Most methods proposed so far in this context generally differ in the way they promote and measure independence. In the remaining of this paper, we will primarily focus on the separation of sparse sources.

A. Sparse blind sources separation

Sparkled by advances in harmonic analysis and applied mathematics, the sparse modeling of signals [START_REF] Elad | Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing[END_REF] has attracted a lot of interest and has proved to be effective at solving a very large range of inverse problems: denoising, deconvolution [START_REF] Starck | Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity[END_REF], compressive sensing reconstruction [START_REF] Candès | People hearing without listening: An introduction to compressive sampling[END_REF], inpainting [START_REF] Fadili | Inpainting and zooming using sparse representations[END_REF] to only name a few. In this context, it is assumed that the each source {s j } j=1,••• ,n can be sparsely decomposed in some basis, waveform dictionary or signal representation Φ:

s j = α j Φ. (2)
The sparsity of s j in the signal representation Φ is customarily described with two different models:

• Exact sparsity model: each source {s j } j=1,••• ,n is assumed to have a small number of nonzero coefficients in Φ. The set of nonzero coefficients, or support of the source j, is denoted by Ω j .

• Approximate sparsity model: natural signals do not generally have an exactly sparse distribution in Φ. They rather exhibit an approximately sparse distribution: most of the entries of the expansion coefficient vectors {α j } j=1,••• ,n are zero or with negligible amplitudes and only a few take significant amplitudes.

Most natural signals exhibit such sparsity property in adequately chosen signal representations Φ. Examples of standard signal representations include the wavelets [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], [START_REF] Starck | The undecimated wavelet decomposition and its reconstruction[END_REF], the curvelets [START_REF] Fadili | Curvelets and ridgelets[END_REF] or even adaptively learned signal representations [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. The sparsity property is particularly appealing since few significant expansion coefficients of s in Φ efficiently or compressively encode its information content. Sparsity has been exploited in recent sparse blind source separation methods (see [15]- [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF]).

In the framework of ICA, sparse and independent sources have been studied in [15], [START_REF] Koldovsky | Efficient variant of algorithm fastica for independent component analysis attaining the cramer-rao lower bound[END_REF]. The RNA method [15] looks for sources by enforcing their sparsity in an orthonormal basis. The Efficient FastICA (EFICA) algorithm described in [START_REF] Koldovsky | Efficient variant of algorithm fastica for independent component analysis attaining the cramer-rao lower bound[END_REF] is an improvement of the acclaimed FastICA algorithm [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF] which is specifically adapted to retrieve sources with generalized Gaussian distribution; this includes sparse sources. The compressibility of the signal allowed by sparse representations has been exploited in the Sparse Component Analysis (SCA) algorithm introduced in [START_REF] Li | Underdetermined blind source separation based on sparse representation[END_REF], [START_REF] Gribonval | A survey of sparse component analysis of blind source separation: principles, perspectives, and new challenges[END_REF], [START_REF] Lesage | Underdetermined source separation: comparison of two approaches based on sparse decompositions[END_REF]. In the framework of SCA, the sources are assumed to have nearly or exactly disjoint supports Ω j : for each source, a few entries are nonzero and each active entry is active in only a single source. Relaxing this assumption has been proposed assuming that there exist entries where only one source is active. Assuming that the sparsity patterns of the sources are fully or partially disjoint has inspired several sparse BSS techniques for the processing of time series and audio signals [START_REF] Jourjine | Blind separation of disjoint orthogonal signals: demixing n sources from 2 mixtures[END_REF]- [START_REF] Abrard | A time-frequency blind signal separation method applicable to underdetermined mixtures of dependent sources[END_REF]. However these assumptions are seldom verified with more complex data, especially in imaging : i) signals or images are generally not exactly but rather approximately sparse in Φ, thus entailing that there is no sample where the sources vanish, ii) there is generally no sample where each source is active alone. These assumptions have been further relaxed to approximately sparse sources in [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF] with the introduction of morphological diversity (MD). This concept rather states that sources with different morphologies are unlikely to have similar large amplitude coefficients in Φ. For instance, sparse and independently distributed sources are very likely to verify the MD assumption. Since high-amplitude entries of the sources in Φ encode the most salient morphological features of the sources, sources which do not share similar morphologies will not therefore activate similar expansion coefficients in Φ. Building on the concept of morphological diversity, Generalized Morphological Component Analysis (GMCA) has been showed to be an effective BSS method [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF]. Recent developments in sparse BSS have mainly emphasize on solving non-negative matrix factorization problems with sparsity constraints on the sources [START_REF] Kim | Group Sparsity in Nonnegative Matrix Factorization[END_REF]- [START_REF]Analysis and Synthesis Sparsity in NMF[END_REF]. In the present paper, we focus on the general setting where the mixing matrix and the sources are not necessarily non-negative.

B. Real-world data are often partially correlated

Dealing with real-life signals reveals that the actual sources are seldom perfectly uncorrelated: neither statistical independence nor morphological diversity are valid assumptions. Most classical source separation methods may not work properly for dependent sources and would likely yield erroneous results.

In the field of spectroscopic data processing, the sources of interest are often partially correlated. In this particular field of research, the development of dedicated source separation techniques has already been considered assuming strong assumptions about the sources to be retrieved [START_REF] Naanaa | Blind source separation of positive and partially correlated data[END_REF], [START_REF] Sun | Postprocessing and sparse blind source separation of positive and partially overlapped data[END_REF]. In [START_REF] Naanaa | Blind source separation of positive and partially correlated data[END_REF], the separation of partially correlated signals is made possible by assuming that for each source there exists at least one non-zero entry that vanishes for the other sources; this is obviously reminiscent of the disjoint source support introduced in the framework of SCA. Together with the non-negativity of the sources and the mixing matrix, the fact that the sources are assumed to have partially disjoint supports dramatically simplifies the separation task. Indeed, in such a setting, each column of the matrix are colinear with some of the columns of the data X. Thanks to the non-negativity assumption, the data are positive linear combinations of positive sources. From a geometrical point of view, the columns of the mixing matrix A are the vertices of the conical hull generated by the columns of X (see [START_REF] Donoho | When does non-negative matrix factorization give a correct decomposition into parts[END_REF], [START_REF] Gillis | The why and how of nonnegative matrix factorization, ser. Machine learning and Pattern recognition series[END_REF] for more precisions). This geometrical property has been exploited to design algorithms that seek the simplex with minimal volume that encloses the data samples [START_REF] Chan | Convex analysis based minimum-volume enclosing simplex algorithm for hyperspectral unmixing[END_REF]. This framework has been showed to be well adapted to the blind separation of very sparse and non-negative sources. This assumption has been somehow relaxed to model sources that are composed of sparse spikes and a non-negative baseline in [START_REF] Sun | Postprocessing and sparse blind source separation of positive and partially overlapped data[END_REF]. In this case, the sources are modeled as a combination of two components: i) a sparse spikes train and ii) a smooth continuum, where the sparse spike component verify the source support assumption made in [START_REF] Naanaa | Blind source separation of positive and partially correlated data[END_REF]. Partially correlated sources are often found in other seemingly different scientific domains such as astrophysics. As an illustration, the left panels of Figure 1 show simulations of two components which originate from distinct astrophysical phenomena. From a mathematical perspective, these two components are characterized by rather different columns of the mixing matrix (in physics, this translates into different electromagnetic spectra), and by somehow different spatial distributions. As a consequence, they can be physically and mathematically considered as sources. These components are generated by different physical mechanisms (e.g. synchrotron emission, Bremsstrahlung, thermal emission, etc.) and are therefore considered as distinct components or sources by astrophysicists. However these signals depend on common physical parameters (e.g. gas density, temperature, etc.) which may make them partially correlated. We refer the interested reader to [START_REF] Bouchet | Foregrounds and cmb experiments: I. semi-analytical estimates of contamination[END_REF], [START_REF] Delabrouille | The pre-launch Planck Sky Model: a model of sky emission at submillimetre to centimetre wavelengths[END_REF] for more details about the astrophysics in the microwave wavelengths. The right panel of Figure 1 reports the values taken by the correlation coefficients between these two astrophysical components in the wavelet domain. This particular example also highlights a second origin for signal partial correlation: chance-correlations. Indeed, in the finite sample case, single realizations of theoretically independent or uncorrelated sources are likely to have non-zero correlation. In Figure 1, the number of available samples decreases at low-frequency, which plays a non-negligible role in the increase of the correlation coefficients at large-scale. In the general setting, such as in imaging science, the limitations of standard methods are the following: Figure 1. Left : samples of partially correlated sources in astronomical imaging: simulations free-free and spinning dust emissions as observed by the WMAP probe. Right : correlation coefficient of astrophysical components per wavelet scale.

• Sparsity: standard methods assume that the sources are sparse and non-negative in the sample domain.

However, it is more generally appropriate to model the sparsity of the sources in some transformed domain Φ, where non-negativity does not hold.

• Non-negativity: the non-negativity of both the sources and the mixing matrix, which is a valid assumption in a wide range of applications, does not hold in some applications (for instance in astrophysics, see [START_REF] Bobin | SZ and CMB reconstruction using Generalized Morphological Component Analysis[END_REF]).

• Sources with partially disjoint supports: assuming that, for each source, there exist at least one data sample where a single source is active is a strong assumption in general. In the general setting, the sources are approximately sparse in a transformed domain. As a consequence, all the sources have non-vanishing entries.

The goal of the present paper is to introduce a novel blind source separation algorithm that relaxes these assumptions in the case of sparse and partially correlated sources.

Contribution

In this article, building upon sparsity and the concept of morphological diversity, we introduce a novel blind source separation technique for the separation of sparse and partially correlated sources from noisy measurements. This new technique coined Adaptive Morphological Component Analysis (AMCA) makes no assumption about the non-negativity of the sources or their mixtures. This makes this method suitable for a very wide range of separation problems especially in imaging. The paper is organized as follows: Section II reviews the basics of the GMCA algorithm, its limitations in the case of partially correlated sources, and how they can be alleviated. Section II-B introduces the AMCA algorithm. Section III reports results of extensive evaluations of performances the AMCA algorithm with respect to standard sparse BSS algorithms. Finally, Section III-D illustrates the performances of AMCA in the context of astrophysical source separation.

II. ADAPTIVE SPARSITY AND SEPARATION OF PARTIALLY CORRELATED SOURCES

A. Morphological diversity and the separation of partially correlated sources

Building on the concept of morphological diversity, the GMCA algorithm (Generalized Morphological Component Analysis -see [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF]) aims at retrieving the sparse sources and the mixing matrix by solving the following optimization problem :

min A,S n j=1 λ j s j Φ T p + 1 2 X -AS 2 F , (3)
where the first term penalizes non-sparse solutions; the quadratic part is a classical data-fidelity term. The Frobenius norm is defined as follows: Y 2 F = Trace YY T . The sparsity level is measured by the p norm of the sources with p ≤ 1. We generally choose either p = 1 or p = 0; we refer the interested reader to [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF] for a detailed discussion of the choice a particular p norm for the sparsity penalty. The 1 norm is particularly appealing since it makes the estimation of S a convex optimization problem when A is fixed. However, as discussed in [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], it tends to produce biased estimates of the sources. In the remaining of this paper, we will therefore choose p = 0 which has been showed to provide the best separation results in [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF]. The problem in Equation ( 3) is classically tackled by using an alternate projected least-square algorithm. The GMCA algorithm processes by iteratively alternating the following steps :

• Step 1 : Estimation of S for fixed A:

min S n j=1 λ j s j Φ T 0 + 1 2 X -AS 2 F .
• Step 2 : Estimation of A for fixed S: min

A X -AS 2 F .
Unless A and Φ are orthogonal matrices, the solution to the problem in Step 1 does not admit a closed-form expression. In the GMCA algorithm, this step is rather approximated with a projected least-square estimator :

∀j = 1, • • • , n; s j = H √ λj [A + X] j Φ T Φ,
where A + = (A T A) -1 A T is the Moore pseudo-inverse of the matrix A and [Y] j stands for the j-th row of the matrix Y. The operator H √ λj is the well-known hard-thresholding operator with threshold λ j . For any vector y, this operator performs on each entry y[t] as follows:

H µ (y[t]) = y[t] if |y[t]| > µ 0 otherwise
The second step of the GMCA algorithm is a more classical least-square problem which admits a closed-form solution :

A = XS T (SS T ) -1 .
a) Sparse BSS and partially correlated sources: In this section, we discuss how sparse and partially correlated (s.p.c.) sources can be estimated using sparse blind source separation. For the sake of clarity, the transform Φ is fixed to the identity matrix. In essence, it has been emphasized in [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF] that, in the low-noise limit, solving the problem in Equation 3 with GMCA is somehow similar to finding the sources which are jointly the sparsest possible. More precisely, assuming that the 1 norm is used to measure the sparsity level of the sources, sparse BSS tends to estimate the sources and the mixing matrix so that the sources are enclosed in the 1 ball with the smallest radius. It is important to notice that, for sparse sources with independently distributed entries, any non-trivial mixture of the sources has, with high probability, a higher 1 norm (see [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF]); it will therefore be enclosed in a 1 ball with larger radius. This observation motivates the design of a separation procedure that looks for the sparsest sources as performed by the GMCA algorithm.

In the case of s.p.c. sources, this separation procedure may no longer be viable. As an illustration, Figure 2 shows an example of two s.p.c. sources. In brief, the nonzero entries of the sources are independently distributed with the exception of a few samples -labeled with red and blue dots -which are shared by the two sources. Figure 3 displays the scatter plot (i.e. this plane features the samples of s 1 in abscissa and s 2 in ordinates) of the two sources estimated by looking for the sparsest sources from two mixtures as is done in the GMCA algorithm. The estimated sources are enclosed in the 1 ball with the smallest radius (dotted purple 1 ball). However they do not match the sought after sources; these would actually be enclosed in the dotted red 1 ball. The black arrow features the axis along which the abscissa and ordinates should be aligned if the right sources were correctly estimated. This simple toy-example highlights the limits of any standard sparse BSS method: the presence of significantly active samples in more than one source entails that the sparsest solution might not be the correct solution.

These empirical observations raise the following crucial question: which are the sources' samples that are the most relevant or informative for separation ? Sparse BSS methods will generally be handicapped by the presence of samples which simultaneously have a significant amplitude in more than one source. For instance, in Figure 3, the sample located at the position (-0.8, 0) takes a large amplitude in both sources. These samples are deemed non-discriminant and therefore of poor relevance for the estimation of the sources and the mixing matrix. In the rest of the paper, we will assume that the samples of the sources can be decomposed into two distinct subsets: Ω c will define the set of the sources' samples where morphological diversity does not hold: these samples are likely to take significant amplitudes in more than one source. The set Ω will denote the complement of Ω c ; morphological diversity is assumed to hold true in this set. Ω will therefore contain the samples that are the most discriminant and, subsequently, the most relevant for the separation. If the subset Ω were known, the sources could be estimated in a straightforward way by restricting the application of any sparse BSS method to Ω . In the framework of GMCA, one would trivially substitute step 2 with the following :

min This assumption is unfortunately unrealistic in practice. The set Ω has therefore to be learned along with the sources and the mixing matrix. It is important to note that attempting to estimate Ω as an extra unknown in Equation ( 3) would lead to a cumbersome combinatorial problem. Since it cannot be reasonably estimated directly as an extra parameter, we rather advocate the use of a more flexible way to privilege discriminant samples, which are likely to belong to Ω . The problem in Equation 4 is equivalent to a weighted least-square problem of the form:

A [X -AS] Ω 2 F , (4)
min A Trace (X -AS)W(X -AS) T , (5)
where W is a diagonal matrix with ones for samples in Ω and zeros otherwise. Since Ω is unknown, we propose relaxing this strategy with an adaptive re-weighting procedure. For that purpose, the mixing matrix will be estimated in a similar manner as in Equation ( 5). Each sample is given a weight which measures their discriminant level. These weights are encoded in the diagonal matrix W. We highlighted previously that standard BSS methods are hindered by samples with large amplitudes in more than one source. Conversely, the most discriminant samples will be the ones that mainly take a large amplitude in a single source. These properties suggest that discriminant samples will be traced by the sparsity level of the columns of the matrix S. Consequently, we define the weight vector w q ∈ R T as follows:

∀t = 1, • • • , T ; w q [t] = 1 [S] t q + , (6)
where [S] t is the t-th column of S. The scalar may be required when there exist columns of S with vanishing q norm; its value is generally small, typically = 10 -6 . The parameter q will be chosen in the range [0, 1] so as to measure the sparsity level of the samples across the sources. Other choices for the samples' weights can be considered. However, this choice is particularly interesting as : i) it accounts for the natural sparsity of the sources, ii) the samples in Ω are indeed sparser across the sources than samples which are common to all sources and iii) this weighting strategy is inversely proportional to the amplitude of the columns of S. The former entails that it penalizes more samples of the sources which are common to several sources and have a high amplitude in the meantime. This type of samples are generally the most detrimental to sparse BSS methods when the sources are partially correlated.

B. AMCA: Adaptive Morphological Component Analysis

In this section, we introduce a novel blind source separation algorithm, which we called AMCA (Adaptive Morphological Component Analysis), to specifically deal with s.p.c. sources, that builds upon the GMCA algorithm. In the previous section, we advocated the use of a weighting strategy to privilege/penalize data samples based on their ability to discriminate between s.p.c. sources. This formally amounts to estimating the mixing matrix A and the sources S by solving the following minimization problem:

min A,S n j=1 λ j s j Φ T p + 1 2 Trace (X -AS)W(X -AS) T , (7)
where the matrix W stands for the weighting matrix. This problem takes a simpler formulation in the transformed domain. For the sake of simplicity, we will first consider that Φ is an orthogonal matrix. We will further define X Φ = XΦ T and S Φ = SΦ T , which respectively stand for expansion coefficients of the data and the sources in the transformed domain. Thanks to this re-parameterization, and assuming that p = 1, the problem in Equation ( 7) can be recast as follows:

min A,SΦ ΛS Φ 1 + 1 2 Trace (X Φ -AS Φ )W Φ (X Φ -AS Φ ) T , (8)
where W Φ = ΦWΦ T , and

Λ = diag(λ 1 , • • • , λ j , • • • , λ T
) is a diagonal matrix that stores the regularization parameters. According to the sparse modeling developed in the previous section, the partial correlation of the sources is modeled as a correlation of certain entries of the sources in the transformed domain. Subsequently, the matrix W Φ is diagonal and it can be defined as W Φ = diag (w q w q ), where the operator is the Hadamard product. From Equation (6), the weight vector is defined as:

∀t = 1, • • • , T ; w q [t] = 1 [S Φ ] t q + . ( 9
)
In the next, we will write denote the weighting matrix as W Φ (q) to highlight its dependency on the parameter q. The optimal weight vector w q clearly depends on the unknown true sources. In practice, one has generally only access to an imperfect estimate ŜΦ ; the weight vector w q and the weight matrix W Φ (q) will be computed from some current estimate ŜΦ in the proposed algorithm. For a fixed weight matrix W Φ , minimizing the problem in Equation 8 can be tackled by alternately updating A and S Φ in the spirit of Block Coordinate Relaxation(BCR -see [START_REF] Bruce | Block coordinate relaxation methods for nonparametric signal de-noising[END_REF], [START_REF] Tseng | Convergence of a Block Coordinate Descent Method for Nondifferentiable Minimization[END_REF]). This procedure then turns to alternately solving sequences of convex minimization problems in place of a globally non-convex problem.

Updating S Φ for a fixed mixing matrix A: Assuming that the mixing matrix A is fixed, estimating the source matrix S Φ is performed by looking for the unique minimizer of the following convex problem:

min SΦ ΛS Φ 1 + 1 2 Trace (X Φ -AS Φ )W Φ (X Φ -AS Φ ) T , (10)
The equation to be minimized can be decomposed into two terms: i) a quadratic and differentiable data fidelity term, and ii) a convex but non-differentiable 1 penalty. Fortunately, the 1 is a simple and proper convex function that admits a proximal operator (we refer the interested reader to [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], [START_REF] Parikh | Proximal algorithms[END_REF] for more details about proximal calculus). Furthermore, the quadratic term is differentiable and its gradient is L-Lipschitz with constant L = AW Φ A T 2 . This entails that the problem in Equation ( 10) can be solved exactly using techniques such as the Forward-Backward splitting algorithm [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. This optimization strategy has been used in [START_REF] Rapin | Sparse and Non-Negative BSS for Noisy Data[END_REF] for solving sparse non-negative matrix factorization problems. However, it has the strong disadvantage of dramatically increasing the computational cost of the algorithm: update S Φ would require resorting to efficient but costly iterative algorithms. Furthermore, each time the source matrix S Φ is updated in the algorithm AMCA algorithm, it is fully re-estimated; therefore, it may not be necessary to update with high precision S Φ at each step of AMCA. Subsequently, we rather resort to a simpler but approximate minimization scheme in the spirit of the GMCA algorithm. More precisely, finding the critical point Equation ( 10) can be done by deriving its subdifferential with respect to S Φ [START_REF] Rockafellar | Convex analysis, ser. Princeton Landmarks in Mathematics and Physics[END_REF] and finding its roots:

0 ∈ Λ ∂ ΛS Φ 1 -A T X Φ W Φ + A T AS Φ W Φ , (11)
where we used the fact that W Φ is diagonal. The term ∂ ΛS Φ 1 stands for the subdifferential of the matrix 1 norm taken at ΛS Φ ; it is well known that the subdifferential of the matrix 1 norm is defined as follows:

∂ Y 1 = Z ∈ R n×T Zj, k = sign(Y j,k ), Y j,k = 0 Zj, k ∈ [-1, 1], Y j,k = 0 .
The root of the subdifferential in Equation ( 11) does not admit an explicit formulation unless A T A is diagonal, which would make Equation [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] separable with respect to each individual source. In such a case, each source S Φj can be estimated independently and takes the following closed-form expression:

∀j = 1, • • • , n; S Φj = S µj 1 a j 2 2 a j T X Φ , (12)
where S µj stands for the soft-thresholding operator (i.e. the so-called proximal operator of the 1 norm) and µ j = λj a j 2 2 WΦ j,j . It is important to note that, in this setting, this is equivalent to:

S Φ = S µ (A T A) -1 A T X Φ , (13)
where µ = [µ 1 , • • • , µ n ] and (A T A) -1 A T X Φ is precisely the least-squares estimate of the sources. In other, this procedure amounts to thresholding the least-square estimate (i.e. the minimum of the quadratic data fidelity term), which thus dramatically reduces the computational burden of this step in the AMCA algorithm. However, this equivalence holds whenever A is orthogonal (i.e. A T A is diagonal), which is not always a valid assumption. As a consequence, updating S Φ according to Equation ( 16) provides a rough proxy in the general case which might prevent the algorithm from convergence. This will be discussed later on in this section.

Updating A for a fixed source matrix S Φ : The mixing matrix A is updated assuming the source matrix S Φ is fixed, which is performed by tackling the following convex problem:

min A Trace (X Φ -AS Φ )W Φ (X Φ -AS Φ ) T . ( 14
)
This is a standard reweighted quadratic problem, which admits the following unique minimizer:

A = X Φ W Φ S Φ T S Φ W Φ S Φ T -1 (15)
The AMCA algorithm: In the spirit of BCR, the AMCA algorithm mainly alternates by updating S Φ and A according to Equations ( 16) and (15). However, the overall estimation problem of Equation ( 7) is nonconvex; this means that this estimation procedure is likely to be prone to trapping in local minima. This will highly depend on the initialization of the AMCA algorithm. In [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF], we showed that thanks to the sparsity of the sources and more specifically to the morphological diversity assumption, starting from large regularization parameters {µ j } j=1,••• ,n (or equivalently large thresholds) and decreasing their values towards a final noise-dependent threshold greatly helped improving the robustness of the algorithm GMCA with respect to local minima. The same procedure is implemented in the AMCA algorithm: the {µ j } j=1,••• ,n are first set up to some large initial values {µ (0) j } j=1,••• ,n and then decreased towards final values {µ

(f ) j } j=1,••• ,n .
As pointed out earlier in this Section, the weight matrix W Φ is computed from some guess for the true sources S Φ . Since one has only access to estimated sources, and these estimates are expected to improve at each step of the algorithm, the weight matrix is updated at each iteration k according to the current estimate of the sources S Φ (k) . The choice of the parameters are discussed in Section II-C. In [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF], different thresholding strategies have been considered, namely soft and hard-thresholding. In contrast to hard-thresholding, the main drawback of the soft-thresholding is that it is prone to biasing the estimated sources. In practice, substituting the soft-thresholding with a hard-thresholding in the GMCA algorithm provided a substantial improvement of the separation performances. As a consequence, the update step described in Equation ( 16) is replaced by :

S Φ = H µ (A T A) -1 A T X Φ , (16)
where H µ stands for the hard-thresholding operator with thresholds µ

= [µ 1 , • • • , µ n ].
The resulting algorithm has been named Adaptive Morphological Component Analysis (AMCA). The adaptivity of AMCA stems from its attempt to estimate a weight matrix W Φ that better favors discriminant entries. In other words, this can be interpreted as adaptively learning which samples are more likely to verify the morphological diversity. The AMCA algorithm is detailed below:

I. Transform the data X Φ = XΦ T II. Initialize A (0) , S (0)
Φ , {µ

j }j=1,••• ,n and q (0) = 1 III. While k < Pmax, (0)
1 -Update S assuming A is fixed :

S Φ = A (k) T A (k) -1 A (k) T X Φ ∀j = 1, • • • , n; [S Φ (k+1) ]j = H µ (k) j (S Φ ),
where [S Φ (k+1) ]j is the t-th column of S Φ (k+1)

2 -Update the weight matrix :

W Φ (k) q = diag   1 [SΦ (k+1) ] t q (k) +   , where S Φ (k+1) t is the t-th column of S Φ (k+1)
3 -Update A assuming S is fixed : On the convergence of the AMCA algorithm: The problem in Equation ( 8) is not convex; therefore, only convergence to a critical point can be expected. The design of the AMCA algorithm has been inspired by Block Coordinate Relaxation [START_REF] Bruce | Block coordinate relaxation methods for nonparametric signal de-noising[END_REF]. In [START_REF] Tseng | Convergence of a Block Coordinate Descent Method for Nondifferentiable Minimization[END_REF], Tseng precisely proved the convergence of BCR for the minimization of nondifferentiable and nonconvex cost functions, with a special application to sparsity-constrained blind source separation problem. According to [START_REF] Tseng | Convergence of a Block Coordinate Descent Method for Nondifferentiable Minimization[END_REF], alternating the minimization steps in Equations ( 10) and ( 14), for fixed parameters {µ j } j=1,••• ,n and W Φ , converges to a critical point of the problem in Equation ( 8). Firstly, decreasing the thresholds {µ j } j=1,••• ,n is a strategy that has been proposed in [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF] to improve the robustness of the GMCA algorithm to spurious local minima. In the field of optimization, this procedure is reminiscent of the fixed point continuation technique [START_REF] Hale | A Fixed-Point Continuation Method for l1 -Regularized Minimization with Applications to Compressed Sensing[END_REF], which has been proposed to speed up the minimization of 1 -penalized least-squares problems. The convergence of the AMCA algorithm would be guaranteed as long as steps ( 10) and ( 14) are alternated until convergence for each value of {µ j } j=1,••• ,n . The thresholds are however updated at each iteration of the AMCA algorithm, which helps speeding up the algorithm but might prevent it from convergence. Secondly, weight matrix W Φ is updated at each iteration, which also might prevent the AMCA algorithm from convergence. Lastly, in the spirit of reweighted 1 techniques [START_REF] Candès | Enhancing Sparsity by Reweighted 1 Minimization[END_REF], the weights W Φ are updated based on the current estimate of the sources S Φ . If this strategy is a well motivated heuristic, the convergence of the AMCA algorithm is not theoretically grounded. In numerical experiments, the AMCA tends to stabilize after a given number of iteration. As an illustration, Figure 4 displays the evolution of the signal-to-distortion ratio (SDE -see Section III for more details) during a run of the AMCA algorithm on simulated sources, which are described in Section III. Our intuition is that, after a given number of iterations, the parameters {µ j } j=1,••• ,n and W Φ are likely to evolve slowly enough to make the AMCA stabilize. In practice, the AMCA algorithm converges in a few hundreds of iterations. Subsequently, choosing the number of iterations P max = 500 provided good results.

A (k+1) = X Φ W Φ (k) q S Φ (k+1) T S Φ (k+1) W Φ (k) q S Φ (
From orthogonal transforms to tight frames: Previously, the derivation of the AMCA algorithm has been carried out assuming that the signal representation Φ is an orthogonal transform; i.e. it verifies: Φ T Φ = ΦΦ T = I. This property allows to perform AMCA in the transformed domain. However, it has long been emphasized that appropriately chosen redundant transforms generally allow for extra degrees of freedom which help producing sparser representations. In the context of sparse BSS, the use of redundant transforms has been shown to yield substantial separation improvements (see [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF]). In the numerical experiments of Section III, we make use of redundant wavelets to sparsely represent the sources. Such type of wavelets have the property of translation-invariance [START_REF] Coifman | Translation-invariant de-noising[END_REF], [START_REF] Starck | Wavelets, ridgelets and curvelets on the sphere[END_REF]. These transforms are non-orthogonal; they are tight frames which verify the exact reconstruction property but for which ΦΦ T = I. In practice, this means that S Φ = SΦ T is a sparse representation of the signal S but not generally the sparsest, if one exists (see [START_REF] Elad | Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing[END_REF]). Fortunately, since the Gram matrix ΦΦ T is diagonally dominant for most tight frames of interest in imaging, a first order approximation amounts to resorting to the same changes of variables S Φ = SΦ T and X Φ = XΦ T . This allows to apply the AMCA algorithm in the transform domain in the tight frame case which greatly lower the computational cost of the algorithm.

C. Choosing the parameters in AMCA

Choice of the thresholds {µ j } j=1,••• ,n : It has been emphasized in [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF] that the choice of the thresholds {µ j } j=1,••• ,n plays a crucial role in the separation process. From the concept of morphological diversity, the high-amplitudes entries of the sources are likely to belong to a single source. In the case of statistically independent sources, this assumption holds true, with high probability, for high enough amplitudes. This makes the high-amplitude entries the most discriminant samples of the sources. The hard-thresholding operator H µj selects entries with amplitudes higher than √ µ j and rejects the others. In the AMCA algorithm, the thresholds {µ j } j=1,••• ,n are first initialized to high values; generally each parameter µ j is set to the maximum amplitude of the initial guess of the j-th source s (0) j . Their values then decrease towards a final threshold which depends on the noise level. In practice, the final thresholds are fixed to 3 σ j where σ j stands for the standard deviation of the noise which contaminates the j-th sources. Whenever the level of the noise is unknown, it is estimated empirically using the Median Absolute Deviation estimator (MAD). Assuming that the noise is additive white Gaussian, this procedure guarantees, with high probability, that no noise sample is selected. Such kind of coarse-to-fine procedure has some interesting properties : i) the sources are first estimated from the most discriminant entries (i.e. the high-amplitude entries of the sources), ii) in the spirit of simulated annealing, the variation of the thresholds makes AMCA less prone to trapping in local minima and iii) the thresholding helps preventing AMCA from incorporating noise entries which makes it more robust to noise.

Choice of the weight parameter q: The AMCA algorithm rely on a re-weighting procedure that penalizes/favors certain entries of the estimated sources. The weights are function of the q norm of the columns of S. They somehow measure the activity of each sample across the sources. Intuitively, choosing a low value for q seems quite natural since it yields more contrast between sparse and non-sparse columns of S. This argument would make perfect sense if the true sources were known. A trade-off has to be made between the two following options : i) large values for q might lead to an under-penalization of less discriminant entries and ii) small values for q provide a larger penalization of non-sparse entries of S, which is desirable to efficiently separate s.p.c. sources. However, at the beginning of the AMCA algorithm, one has only access to imperfect, if not erroneous, estimates of S. In this case, small values of q might mis-penalize/mis-favor entries of S which can eventually hamper the separation process. Alleviating this dilemma is made by starting with a high value for q -typically 1 -and then decreasing it, at each iteration, towards some final value q f . Several values for q f have been tested; it turns out that choosing q f = 0.01 leads to a good trade-off for all the experiments we carried out. Smaller values for q f did not bring any noticeable improvement.

III. NUMERICAL EXPERIMENTS

A. Experimental protocol

In this section, we evaluate the performances of the AMCA algorithm with respect to standard sparse BSS methods : i) the GMCA algorithm [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], ii) RNA (Relative Newton Algorithm) [15] and iii) EFICA (Efficient FastICA) [START_REF] Koldovsky | Efficient variant of algorithm fastica for independent component analysis attaining the cramer-rao lower bound[END_REF]. It is important to recall that, to our knowledge, the blind separation of s.p.c. sources has never been dealt with in the general case where: i) neither the mixing matrix nor the sources are assumed to be non-negative, ii) the sources do not have (quasi)-disjoint supports. Relaxing these assumptions is interesting so as to tackle the separation of complex data which are approximately sparse in general signal representations Φ (e.g. Fourier, wavelets, curvelets, etc.) but where non-negativity is certainly not a valid assumption.

In this first study, our prime intention is to study the performances of standard sparse BSS methods and the AMCA algorithm to estimate s.p.c. sources. For that purpose, we will first make use of synthetic but rather complex data; this makes Monte-Carlo simulations possible in various settings. To mimic spectroscopy-like data, the (s.p.c.) sources are generated according to the so-called (SPC) model introduced below :

• The sources {s j } j=1,••• ,n are K-sparse signals. This means that each source s j has only K non-zero entries out of T entries. The activation of the entries of the sources will be distributed according to a Bernoulli process π with parameter 0 ≤ ρ ≤ 1 so that :

π = 1 with probability ρ, 0 with probability 1 -ρ. ( 17
)
• For each source, L ≤ K active entries are shared by all the sources and the remaining K -L are independently drawn. Let us recall that Ω c is the set of the L entries which are common to all the sources; the distribution of these entries is Gaussian with mean 0 and variance τ 2 . In this model, common active entries do not have correlated amplitudes.

• The K -L nonzero entries which belong to Ω j \Ω c are are independently distributed according to a Gaussian law with mean 0 and variance 1. The parameter τ defines some dynamic range between the amplitude of common and independent samples of the sources.

• the sources are convolved with a Laplacian kernel of full width at half height (FWHM) equal to 15. An example of 2 simulated sources (resp. 2 mixtures) is shown in the left (resp. right) panel of Figure 5. The resulting sources can be sparsely represented in a frame of redundant (translation invariant) 1D wavelets 1 . Each source is made of T = 4096 samples. The matrix will be picked at random from a Gaussian distribution. In the following experiments, the sparsity level K/T is fixed to 0.02. This means that for T = 4096 -which will be the actual value in the next -about 81 spikes will be active in each source.

b) Evaluation criteria: Following [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF], the estimated sources can be decomposed into the combination of contributions of different origins:

s est = s target + s interf + s noise + s artefacts ,
with the following interpretations for the terms:

• s target is the projection of s est on the target (ground-truth) source. In other words, it is the one part of this decomposition which corresponds to what needs to be recovered. The other ones are residues. • s interf accounts for interferences due to other sources.

• s noise is the part of the reconstruction which is due to noise.

• s artefacts stands for the remaining artefacts which are neither due to interferences nor noise.

From this decomposition, one can derive a global criterion named Source to Distortion Ratio (SDR) which contains information about each of the different terms of the decomposition: SDR(s est ) = 10 log 10 s target 2 2

s interf + s noise + s artefacts 2 2 . ( 18
)
In [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], we advocated the use of a global performance criterion that only depends on the mixing matrix. Such a criterion is particularly interesting in the noisy setting as it is not dependent on the denoising effect of the methods to be compared: methods like GMCA or AMCA have the ability to provide denoised estimates of the sources while others are not. For that purpose, the following mixing matrix criterion C A is introduced:

C A = PA + A -I 1 . ( 19
)
where PA + is the pseudo-inverse of the estimated mixing matrix corrected -through P -for the scale/permutation indeterminacies; A is the true mixing matrix. Low values of C A then indicate better separation performances.

In the rest of the paper, all criteria will be computed as the median value of independent Monte-Carlo simulations. Unless specified, the number of simulations for each of the figures was set to 100.

B. Comparisons with standard sparse BSS methods 1) An example of recovered sources: To better visualize how sparse BSS algorithms perform on s.p.c. sources, Figure 6 displays 2 sources out of 10 estimated from 10 mixtures with c = 0.6 and τ = 4. The number of observations m as well as the number of sources are set to 10. The number of samples per source is fixed to 4096. Since c = 0.6, this experiment turns out to be an arduous problem as 60% of the active spikes are shared by all the sources. Figure 6 shows that standard BSS methods all seem to recover some of the most prominent spikes of the sought after source. However they all exhibit a certain level of interferences which indicates poor separation quality. The source estimated by AMCA does not visually exhibit such interferences and perfectly matches the source to be retrieved. 2) Coherence level: In this paragraph, the performances of the AMCA algorithm are evaluated when the level of correlation or more generally coherence, defined by the ratio c = L/K, varies between 0 (the no-correlation case) and 1 (the sources all share the same entries). Its values vary in the range [0, 1]; when c = 0, all the sources are independently distributed. The opposite extreme case c = 1 corresponds to sources exactly the same common active entries which means, according to our model, that the sources are fully correlated sources. The performances of standard sparse BSS methods should decrease for high valued of c. In this experiment, the number of sources and the number of channels are both fixed to 10. The standard deviation of the entries in Ω c is fixed to τ = 4. The noise level is fixed to 120dB. The evolution of the SDR as a function of the coherence level c is displayed in the left panel of Figure 7. First, it is interesting to notice that the performances of the standard BSS algorithms (GMCA, EFICA and RNA) are quickly hampered by the partial correlation between the sources. The GMCA algorithm seems to be less sensitive to the correlation of the sources for c < 0.2. When more than 20% of the entries are shared by all the sources, these algorithms behave equally badly. The AMCA algorithm already shows a much higher SDR even when only a few entries are correlated. When less than 80% of the entries are correlated, the SDR of the AMCA algorithm is higher than 60 dB. Its performances start degrading for c > 0.8 to eventually behave similarly to the standard BSS algorithm when c tends towards 1. The evolution of the mixing matrix criterion is shown in the right panel of Figure 7. One can notice that the mixing matrix criterion with the AMCA algorithm is consistently more than 1 order of magnitudes lower than for the standard BSS methods. Figure 8 displays the average over 100 random simulations of the minimum SDR over the estimated sources; in other words this quantity measures the SDR the most badly estimated source for each simulation. This figure shows that the GMCA algorithm is very good at recovering some of the sources but not all with the same quality. From that respect, the standard BSS methods behave quite similarly : for c > 0.2 they are not able to identify all the sources (the minimum SDR is close to 0 dB). The minimum SDR is 50 dB higher for c < 0.8 with the AMCA algorithm. Consequently, the AMCA algorithm is good at recovering all the sources even when the number of common spikes in the sources is equal to 80%.

3) Dynamic range: A second parameter which is very likely to hamper the performances of standard BSS methods is the dynamic range between between the amplitudes of correlated and non-correlated entries. More precisely, standard sparse BSS algorithms are highly sensitive to the coefficients of the sources which have the highest amplitudes. This makes perfect sense when morphological diversity holds true: the most significant entries of the sources are most likely the most discriminant. Furthermore, we claimed earlier in this paper that the weighting strategy of the AMCA algorithm penalizes more high amplitude coefficients which are shared by several sources. The AMCA algorithm should therefore be less sensitive to the dynamic range between correlated and non-correlated samples of the sources than the standard BSS methods. In this experiment, the number of sources and the number of channels are both fixed to 10. The proportion of correlated entries between the sources is fixed to c = L/K = 0.2. The noise level is fixed to 120dB. Figure 9 features the evolution of the SDR (resp. mixing matrix criterion C A ) in the left (resp. right) panel as a function of the standard deviation of the correlated entries of the sources τ in the interval [0.1, 32]. For τ < 2 (i.e. on average, correlated entries have twice the amplitude of the non-correlated entries), standard BSS methods all behave quite well with SDR values higher than 40 dB. For higher values of τ their performances dramatically decrease towards low SDR values. As expected, the AMCA algorithm is only slightly impacted by the variation of the dynamic range. Even for τ = 32, the SDR of the sources with the AMCA algorithm is higher than 60dB and the corresponding mixing matrix criterion is about 3 orders of magnitude lower than with the standard BSS algorithms.

4) Number of sources:

Recovering a large number of sources is generally a complicated task. Indeed, for a fixed number of samples T , the ability of source separation techniques to disentangle a growing number of sources is rapidly limited. This is true in the case of s.p.c. sources as the number of entries of the sources which are most relevant for the separation are also limited. In the next experiment, the number of sources and the number of channels are both fixed to the same value. The proportion of correlated entries is fixed to c = L/K = 0.2. The noise level is fixed to 120dB and τ = 4. Figure 10 displays the evolution of the SDR on the left panel and C A on the right panel as a function of the number sources in a range 2 ≤ n s ≤ 128. As expected the performances of the standard sparse BSS algorithms, as measured by the SDR, decay as the number of sources to be recovered grows. The AMCA algorithm provides good results (the SDR is higher than 60dB) for less than roughly 64 sources. For a higher number of sources, the performances of the AMCA decrease quite rapidly. However, one has to keep in mind that the number of samples per source is kept fixed to T = 4096 from which 2% are active, that is to say on average 81. Amongst these active samples, only 80% have been drawn independently. This means that each source has about 65 likely discriminant samples out of 4069. When the number of sources increases, the probability for two independently drawn samples to take simultaneously significant amplitudes grows. This might become likely when the number of sources is large which explains why the performances of AMCA decay for n > 64. The evolution of the mixing matrix criterion has to be analyzed with care as the performances of all the methods seem to improve as the number of sources increases while the evolution of the SDR shows a quite different behavior. It has to be noticed that the mixing matrix criterion averages the scalar product of the rows of the inverse of the estimated mixing matrix and the true mixing matrix. In other words, it looks at a quantity that only involves interferences between two sources. On the other hand, the SDR is a global measure of distortion that impacts each of the sources. Precisely, for a fixed value of the mixing matrix criterion, the distortion will increase with the number of sources which explains why the actual SDR of the standard sparse BSS methods decrease with n. C. The AMCA algorithm and noise Discussion about the impact of noise: In this section, we discuss the impact of the re-weighting scheme on the performances of the AMCA in the noisy setting. First, in the AMCA algorithm, the weights w q are estimated from the estimated sources. These sources are obtained via Step 1 of the AMCA algorithm. In the low noise limit, one interesting feature of the proposed re-weighting scheme is that it is inversely proportional to the amplitude of the columns of S. More precisely, this entails that large entries of S which are shared by several sources are more penalized than small entries with the same relative distribution across the sources. Strongly penalizing large and correlated entries is desirable since they are detrimental to the estimation of the mixing matrix and the sources. In the noisy setting, the situation turns out to be rather different since small-amplitude samples are more likely perturbed by noise than large amplitude entries. On the one hand, the proposed re-weighting procedure might be disastrous for the separation of the sources whether they are partially correlated or not. Indeed, since the weights w q are inversely proportional to the amplitude of the columns of the sources, the proposed procedure will tend to favor small entries which are more affected by the presence of noise. On the other hand, Step 1 of the AMCA algorithm rejects entries with amplitudes smaller than some prescribed noise-dependent level. This should lower the impact of noise on the performances of the AMCA algorithm.

Comparisons: In [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], the authors demonstrated that the GMCA algorithm is robust to additive noise contamination. This is especially true whenever morphological diversity holds; in that case the most discriminant sources are the entries of the sources with the most significant amplitudes. It turns out these entries are also the least contaminated by additive noise. In the case of s.p.c. sources, the most discriminant sources are not necessarily the large-amplitude samples. A first consequence is that noise will very likely have a strong impact on the quality of the separation. In this experiment, the number of sources and the number of channels are both fixed to 10. The proportion of correlated entries of the sources is fixed to c = L/K = 0.2, τ is set to 4. The noise level varies in the range [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF]120]dB. The evolution of the mixing matrix criterion is featured in Figure 11. It first reveals that the GMCAbased methods have slightly better performances when the noise level increases. When the noise level tends towards 20dB all these methods tend to all perform badly. More surprisingly, the performances of the AMCA algorithm do not seem to be more impacted by the presence of noise than standard BSS methods. To better understand this behavior, let us recall that, in the AMCA algorithm, the mixing matrix is updated from the current estimate of the sources which is naturally denoised in Step 1 of the algorithm. This means that coefficients of the sources which are of the order of or below the noise level are naturally rejected in Step 1; the estimation of the mixing matrix won't therefore account for these entries. In other words the re-weighting procedure won't "amplify" the impact of the noise. Eventually, for high noise levels, only high amplitude coefficients will be detected. In this regime, the GMCA and AMCA algorithms will be similarly handicapped by correlated large-amplitude coefficients.

D. Application in Astrophysics

Multi-wavelength instruments and sensors have long been used in Astrophysics. However, it is only recently that the design of modern and sophisticated data analysis methods, such as blind source separation, play a significant role in this field. The recent results of the ESA/Planck space mission in 2013 has shed light on the crucial role that blind source separation has played to provide a very accurate map of the Cosmological Microwave Background (CMB) from multi-wavelength microwave data [START_REF] Ade | Planck 2013 results. XII. Component separation[END_REF], [START_REF] Bobin | Joint planck and wmap cmb map reconstruction[END_REF]. In this context, each observation measures the mixture of several astrophysical components including the CMB and galactic foreground emissions to only name two. Galactic emissions can further be decomposed into several physically relevant sources with different spectral contributions. As emphasized in the introduction, these galactic sources are emitted from similar regions of the sky which entails that these components are, by nature, partially correlated. Disentangling between the different galactic sources is a challenging task which has been rarely tackled using blind source separation so far. In the rest of this section, we focus on the estimation of the various galactic sources from CMB-free data -the CMB is customarily removed from the data before studying galactic components. In the range of the microwave wavelengths observed by instruments such as WMAP [START_REF] Bobin | Wmap nine-year cmb estimation using sparsity[END_REF], [START_REF] Bennett | Nine-year wilkinson microwave anisotropy probe (wmap) observations: Final maps and results[END_REF] and Planck, the most prominent galactic emissions are : the free-free emission, galactic synchrotron emission, spinning dust and thermal dust emissions -see [START_REF] Delabrouille | The pre-launch Planck Sky Model: a model of sky emission at submillimetre to centimetre wavelengths[END_REF] for more details about astrophysical emissions in the range of microwave wavelengths and their simulations. In this section, our goal is to study the impact of the partial correlation of astrophysical components on the performances of blind source separation methods. This study will therefore be carried out on simulations of the 5 WMAP [START_REF] Bennett | Nine-year wilkinson microwave anisotropy probe (wmap) observations: Final maps and results[END_REF] observations in the range 23 -94 Ghz generated by the PSM (Planck sky model2 -see [START_REF] Delabrouille | The pre-launch Planck Sky Model: a model of sky emission at submillimetre to centimetre wavelengths[END_REF]). Simulations of these four major galactic emissions at the vicinity of the galactic center are displayed in Figure 12. In these simulations, for the sake of simplicity, the linear mixture model (1) holds true. As displayed in Figure 12, the most prominent emissive areas of all four sources are located in similar regions of the sky. This makes them partially correlated, especially at large scale or low frequencies. The 5 simulated WMAP observations are featured in Figure 13. The sources displayed in Figure 12 have an approximately sparse distribution in a bidimensional translation-invariant wavelet transform [START_REF] Starck | The undecimated wavelet decomposition and its reconstruction[END_REF]; this will be our choice for the signal representation Φ in this section. All the sparse BSS algorithms used in the forthcoming experiments have been performed in the wavelet domain using default parameters. In the context of astrophysics, the impact of noise is usually deemed an important contribution that can hamper the quality of separation. The performances of sparse BSS are therefore evaluated for different levels of noise using the SDR and the mixing matrix criterion C A defined above. The left column of Figure 14 (resp. 15 and 16) features the sources estimated by the EFICA algorithm (resp. RNA and GMCA) for a normalized noise level of 1. The normalized noise level corresponds to the nominal level of the WMAP data. It corresponds to a signal-to-noise ratio of 62dB. This value might seem very large but only large-scale features are emerging above the noise level for mild values of the SNR. With the exception of the estimated free-free emission which shows some spurious point sources residuals that originally belong to the thermal dust emission, all the four sources seem to be correctly estimated. The column on the left of Figure 14 (resp. 15 and 16) features the residual error defined by the difference between the estimated and original sources. Figure 17 displays exactly the same pictures for the AMCA algorithm. The display levels of the residual errors are kept the same to better reveal the differences between the results from all tested algorithms. The level of estimation error provided by the AMCA algorithm seem to be significantly lower for the free-free, spinning dust and thermal dust emissions while tending to be slighter better for the synchrotron emission. Figure 18 shows the evolution of the SDR for the four tested sparse BSS methods when the normalized noise level varies from 0.01 (i.e. the noise level is actually 100 lower than for that real WMAP data) to 2. Each entry of this figure has been obtained as the median value of 100 independent Monte-Carlo simulations. This figure shows that the performances of the standard sparse BSS methods seem to be quite similar at all levels of noise with values of SDR in the range 22 -28 dB. The AMCA algorithm provides estimates with a level of SDR that is almost 10 dB higher in the low noise regime. We already emphasized that the evolution of the SDR yields a partial view of the performances of these algorithms when the noise level increases. In this regime, the mixing matrix criterion C A is better suited to evaluate the quality of separation independently of the level of noise that contaminates the estimates of the sources. The evolution of C A is provided in Figure 18. This figure confirms that the standard sparse BSS algorithms perform in a rather similar manner with values of the mixing matrix criterion of about 0.03. In the meantime, the AMCA algorithm performs slightly better, especially in the low noise regime, with a value of C A that is about 7 times lower in this case. In this high noise level regime, only the large-scale features of the sources are detected. Since these features are likely to be at the origin of the partial correlation of the sources, all algorithms, including AMCA, are handicapped in this setting.

SOFTWARE

Following the philosophy of reproducible research [START_REF] Buckheit | WaveLab and Reproducible Research[END_REF], the codes and algorithms introduced in this article will be available at http://www.cosmostat.org/GMCALab.

IV. CONCLUSION

The blind separation of sparse and partially correlated (s.p.c.) sources is a challenging task which standard sparse BSS methods generally hardly tackle with success. In this article, we introduce a novel sparse blind source separation coined Adaptive Morphological Component Analysis which is designed to retrieve sparse and partially correlated sources. We emphasize that recovering s.p.c. sources can be approached by identifying and privileging the most discriminant entries of the sources. We then propose learning adaptively these entries via a re-weighting scheme jointly with the estimation of the sources and the mixing matrix yielding the proposed AMCA (Adaptive Morphological Component Analysis) algorithm. We further evaluate the performances of the AMCA algorithm with respect to standard sparse BSS methods in various experimental scenarios using Monte-Carlo simulations of s.p.c. sources. The AMCA algorithm is shown to be much less sensitive to partial correlations of the sources even when the sources are highly correlated (i.e. when about 80% of the active samples are common to all the sources). The proposed algorithm is slightly impacted by the dynamic of the amplitudes of correlated and independent entries. It is competitive when the number of sources to be retrieved is large. Finally, we applied the AMCA algorithm to the separation of astrophysical components from simulated microwave data. This application illustrates that the AMCA algorithm is well suited to estimate physical components which are, by nature, partially correlated.

ACKNOWLEDGMENT

This work was supported by the French National Agency for Research (ANR) 11-ASTR-034-02-MultID.

Figure 2 .

 2 Figure 2. A realization of 2 partially correlated sources displayed in blue and green. Dots indicate the active entries of the sources which are common to the sources.

Figure 3 .

 3 Figure 3. Scatter plot of the sources with the sparsest distribution. The boundaries of the 1 ball of the true sources (resp. with the lowest 1 norm) is displayed in purple (resp. in ref). The black arrows show the actual directions along which the scatter plot of the true sources should be aligned.

k+1) T - 1 4-

 1 Decrease each {µ(k) j }j=1,••• ,n and the parameter q IV. Reconstruct the sources S = S Φ Φ.

Figure 4 .

 4 Figure 4. Evolution of the SDR of simulated sources during a run of the AMCA algorithm. In abscissa: iteration number. In ordinate: valued of the SDR.

Figure 5 .

 5 Figure 5. Left: 2 sources generated from the s.p.c. sources model. Right: 2 mixtures of s.p.c. sources.

Figure 6 .

 6 Figure 6. Estimated s.p.c. source. Top left : EFICA, right: RNA. Bottom left: GMCA, right: AMCA.

Figure 7 .

 7 Figure 7. Evolution of the mean SDR (top) and the mixing matrix criterion (bottom) as a function of the coherence level. Each sample is the average over 100 Monte-Carlo simulations.

Figure 8 .

 8 Figure 8. Evolution of the minimal value of the SDR across the sources as a function of the coherence level. sample is the average over 100 Monte-Carlo simulations.

Figure 9 .

 9 Figure 9. Evolution of the mean SDR (top) and the mixing matrix criterion (bottom) as a function of the dynamic range. Each sample is the average over 100 Monte-Carlo simulations.

Figure 10 .

 10 Figure 10. Evolution of the mean SDR (top) and the mixing matrix criterion (bottom) as a function of the number of sources. Each sample is the average over 100 Monte-Carlo simulations.

Figure 11 .

 11 Figure 11. Evolution of CA as a function of the noise level (in dB). Each sample is the average over 100 Monte-Carlo simulations.

Figure 12 .

 12 Figure 12. WMAP data : the 4 input sources.

Figure 13 .

 13 Figure 13. WMAP data : the 5 observations.

Figure 14 .

 14 Figure 14. Left : estimated EFICA sources. Right: estimation error defined by the difference between the estimated and original sources. From top to bottom: free-free, spinning dust, synchrotron and thermal dust.

Figure 17 .

 17 Figure 17. Left : estimated AMCA sources. Right: estimation error defined by the difference between the estimated and original sources. From top to bottom: free-free, spinning dust, synchrotron and thermal dust.

Figure 18 .

 18 Figure 18. Evolution of the mean SDR (top) and the mixing matrix criterion (bottom) as a function of the noise level. Each sample is the average over 100 Monte-Carlo simulations.

For that purpose we used the RICE wavelet toolbox : http://dsp.rice.edu/software/rice-wavelet-toolbox

http://www.apc.univ-paris7.fr/ delabrou/PSM/psm.html


Jason Lambert
email: jason.lambert@cea.fr

Lionel Lacassagne

Gilles Rougeron

Stéphane Le Berre

Sylvain Chatillon

High performances simulation of ultrasonic fields for Non Destructive Testing

Keywords: non-destructive testing, ultrasonic simulation, general purpose processors, graphic processing units, massive parallelism

Ultrasonic imaging is a commonly used method to detect and identify defects in a mechanical part in nuclear applications. Nowadays massively parallel architectures enable the simulation of ultrasonic field emitted by a phased array transducer inspecting a part across a coupling medium. In this paper, regular field computation model will be discussed along its implementations on General Purpose Processors (GPP) and Graphic Processing Units (GPU).

I. Introduction

Inspection simulation is used in a lot of non-destructive evaluation (NDE) application: from designing new inspection methods and probes to qualifying methods and demonstrating performances through virtual testing while developing methods. The CIVA software, developed by CEA-LIST and partners is a multi-technique platform (Ultrasonic Testing (UT), Computed Tomography and Radiographic Testing (CT-RT) Eddy-current Testing (ET), Guided Wave (GW)) used to both analyze acquisitions and to run simulations validated against international benchmark.

In particular, the simulation of ultrasonic field radiated in specimen is widely used in order to design or evaluate probe potential efficiency for a given control. Thus, ultrasonic beam main characteristics such as focal spot or local direction can easily be determined. A common nuclear application is the inspection of steel pipes and nozzles, as illustrated in Figure 1. It shows the simulated ultrasonic fields in a weld between two pipes with a contact probe as opposed to a flexible phased array probe. In the second case the beam is radiated with better directionality and increased amplitude. However, due to the potential complexity of the configurations, current semi-analytical models implemented in CIVA software, based on asymptotic developments, takes minutes to hours to run a simulation. Significant efforts are made to reduce computation times by working both on models and computational aspects. The works reported here ties in with this general approach.

Nowadays, new massively parallel architectures can empower computational software, at the expense of adapting algorithms to the specificities of those architectures to highlight and improve parallel computation steps. The exploitation of both multicore-GPP and GPU capabilities results in a new intensively parallel algorithm of simulation of UT beam. Those architectures have already been used to provide consistent speedup over wave propagation [START_REF] Molero | Accelerating numerical modeling of wave propagation through 2-D anisotropic materials using OpenCL[END_REF], or field modeling for a probe alone [START_REF] Romero-Laorden | Field modeling acceleration on ultrasonic systems using graphic hardware[END_REF] but they have not been applied to UT field simulation yet.

The new model relies on analytical solution to the beam propagation. Its implementations on both architectures use high performances signal processing libraries (Intel MKL on GPP and NVidia cuFFT on GPU). The new model is, however, limited to canonical configurations due to the strict parallelism requirements and to the lack of genericity of analytical beam propagation.

II. Beam propagation modelization for regularity

The ultrasonic field computation relies on the pencil method, a generic approach for heterogeneous and anisotropic structures. By evaluating the ray path of the beam, from the transducer to the observation point, it is possible to evaluate the time of flight and the amplitude of the contribution of the beam using energy conservation principle on the tube. As seen in Figure 2, the beam may propagate through different materials and cross multiple interfaces: its contributions  is determined by the propagation matrix obtained by multiplying the elementary contributions of each section of the pencil with the initial contribution , as shown in equation 1 [START_REF] Gengembre | Pencil method for ultrasonic beam computation[END_REF].

(1)

Figure 2 -Ray tube visualization

In general case, there is no simple solution to determine the ray path from the source point to the computation point. However, in isotropic and homogeneous structures, analytical methods can be used to determine ray path. In direct mode, for standard geometrical surfaces, it is possible to determine a polynomial modelizing the path following Snell-Descartes. The roots of this polynomial, whose degree vary from 4th for planar surfaces to 16th for torical surfaces, correspond to the possible solutions for the ray path. They are determined numerically, through Newton's and Laguerre's Method [START_REF] Lambert | Performance evaluation of total focusing method on GPP and GPU[END_REF]. In the case of half-skip mode, ray path can only be determined analytically for planar surfaces and planar backwalls as illustrated in Figure 3, through two dimensional Newton's method solving.

R1 c1 R2 c2 θ1 θ2 R1 c1 R2 c2 R3 c3 θ1 θ2 θ3 θ2fd

Figure 3 -Ray path determination in direct and half-skip mode

Those methods rely on a numerical resolution of analytical formulas which allow for a greater regularity benefiting to the requirements of massively parallel architectures. Moreover, as this model is dedicated only to a specific set of geometries, the propagation matrix can be fully determined preemptively thus avoiding costly matrix multiplication, for the computed modes. According to the simulated propagation mode (longitudinal or transversal), transmission and reflection Fresnel's coefficients are computed with analytical, specialized, equations.

In each point field, once pencils are computed on the whole transducer, their elementary contributions are summed up to obtain the impulse response of the elastodynamic displacement, as presented in Figure 4.

II. Implementations

In this section, choices of implementation, both on GPP and GPU architectures, will be discussed under the specificities of the hardware. It is noteworthy that the current implementations are restricted to a subset of configurations as this preliminary work is aimed at validating the model and at providing first benchmarking results. The scope of those implementations is planar surfaces and direct mode simulation.

GPP implementation

The GPP implementation is dedicated to take advantage of modern day GPP, composed of multiple general purpose core aimed at executing independent, heavyweight tasks. GPP disposes of two parallelism levels.

• A fine grained parallelism relying on specific SIMD instructions (Single Instruction Multiple Data) which execute the same operation on short vectors (128 to 512 bits). For example, with 128-bit vectors, a SIMD instruction can perform simultaneously four additions on four single-precision floating point numbers (32-bit).

• A coarse grained parallelism relying on multithreading to enable multiple logical tasks to reside simultaneously on the GPP. The OpenMP API is aimed at shared-memory parallelism: it creates a thread per GPP core, each with its own stack where local variables are located but they can also communicate through some shared variable. Its work distribution relies on a succession of sequential sections (with only one active thread) and parallel sections (with all threads active) assembled in a fork-join fashion.

In this implementation, computations are aggregated by coarse step, regrouped over the whole set of points, to even the computation load on the GPP between the cores and in order to maximize the reutilization of the data describing the simulation.

It is noteworthy that both displacement convolution and maximum extraction (relying on the determination of the maximum of the envelope of the displacement modulus) rely on Fast Fourier Transform (FFT). As this operation has been extensively studied and is not the main subject of this work, to benefit from heavily optimized FFT operation, this implementation will use the Intel® Math Kernel Library (MKL). This library offer highly vectorized FFT function to benefits from the SIMD capability of GPP cores.

However, to obtain the best performances, the size of the signal over which the FFT will be applied need to be known in advance: some coefficients can therefore be computed in advance, and then be reused several times over different signals. It is the same notion of plan as that used by cuFFT and other dedicated libraries and it is so-called descriptor in this implementation..

GPU implementation

Due to the high need for computational regularity on GPU, the control flow of this implementation will aim at executing the same task over multiple data and proceed in steps, each corresponding to dedicated kernel. To benefit from highest GPU performances, kernels were developed to use single-precision floating point operations (IEEE-754). To perform signal processing operations, this implementation relies on the optimized cuFFT library which is optimized to perform efficiently on large batches of signals. This library establishes a "plan", consisting in the precomputed coefficients which are applied to the signals. It is most efficient when performed repeatedly on a batch of signals of the same dimensions (number of signals and signal length), to reuse a previously computed plan. The following Figure 6 illustrates the overall required computations steps detailed as follow.

The first step of the computation determines the temporal width of the resulting signal corresponding to the summation of the different elementary contributions obtained from the pencils. It combines two successive kernel calls:

• The first kernel is specialized toward direct mode; threads from a single group work collaboratively to computes the elementary displacements from pencils data before reducing temporal information. Those threads groups are spread across the different points where field is to be computed.

• The second kernel aims at reducing all those temporal information on the whole zone to obtain the global size.

Due to the numerous points where signal should be determined and due to the required size of the signal, it is not possible to compute the whole UT field at once: work is divided on several slices to fit the available memory on a high end consumer grade GPU (1.5GB). The following steps will be repeated on each slice in a loop controlled by the host system.

Loop over points P controlled by the host 3x, on x,y,z Once the maxima are extracted, the host then proceeds with the loop over the next slice of field points.

III. Model validation

In this section a set of configurations, illustrating the simulation capabilities, will be described. The validation will first focus on the impact of the discretization of the probe for pencil construction. Then, once the discretization step is fixed, the validation will address simulation validity. The reference to compare simulations is CIVA 11.0 software. The validated implementations rely on single precision floating point (IEEE-754) for their computations.

Reference configurations description

A set of configurations will be studied, both for model validation and performances. Those consist of planar part made of homogeneous isotropic steel inspected with a linear phased array probe which varies in number of 1x10mm elements. The computation zone consists in 101x101 points, covering a 50x50mm area; it may be in the inspection plane or perpendicular to the beam. The delay law lead to the focusing of the beam in the center of the region. Table 1 summarizes the specificities of each configuration. Similarly, with more samples per elements, the Figure 10 indicates that reducing the discretization width do not reduce the gap between CIVA and this new model. Therefore, in the following paragraphs, 5 samples per elements of the sensor will be used, corresponding to a 1x2mm sampling.

Table 1 -Reference configurations

Model validation

Once that sensor discretization is fixed, the results obtained with this new model can be characterized. On the different datasets, three parameters will be studied: the maximum of amplitude value and the width and height of the focal spot. Those will be compared between the model and CIVA 11.0 reference. The focal spot dimensions are obtained by measuring the distance between the maximum of amplitude and a decrease at -3dB. 2 presents the measurement results for all the studied configurations and the relative error compared to CIVA. First, it is noteworthy that this model slightly overestimates the amplitude obtained by CIVA 11.0; however, in all cases, the discrepancy is less than 0.4 dB (~5%) and remains below the requirements for passing benchmarks against experimental results. Moreover, the focal spot is of the same dimension with both tools, the relative observed error is inferior to 2%. Lastly, results indicate the absence of temporal shift between CIVA and this model (as observed on previous Figures). Thus

Table 2 -Model to CIVA comparison

IV. Model benchmarking

The two implementations (GPU and GPP) have been benchmarked on high-end hardware to aim at a fair evaluation of performances.

GPP benchmarking

Two GPP configurations have been studied: -GPP1 -2x GPP Intel Xeon 5590@3.47Ghz + 24GB of RAM : which disposes of 2x6 hardware core (2x12 logical cores with HyperThreading); -GPP2 2x GPP Intel Xeon E5-2650@2.00Ghz + 24GB of RAM : which disposes of 2x8 hardware core (2x16 logical cores with HyperThreading); Over those two GPP configurations, several tools have been studied through two versions of the implementation: one using a legacy FFT implementation compiled with Microsoft Visual C++ 20101 versus one using the Intel® MKL built with the ICC compiler 2 .

In this paragraph, the study will focus successively, on the scaling capability of this implementation, then on the impact of the hyperthreading before focusing on the impact of the Intel tools.

Figure 13 -Scaling of the application

In those graphs, execution of a L0_64E configuration, with 200x200 field point are studied on GPP1 and GPP2 (16threads). The scaling is then computed by comparing the monothread execution time to the actual time with N threads. This scaling is then presented in the following graph, normalized by the number of threads (a theorical perfect scaling is 1). Compiler : MSVC. The first benchmark aims at measuring the scaling of the application as illustrated in Figure 13. It is noteworthy that the attained speedup is decreasing as the number of threads increases; however the overall speedup is 87% of the maximum. This shows clearly that the elementary impulse response contributions benefits from a great speedup, with a sustained one over 95%. However, the scaling of the signal formation drop to 85% as the number of threads increases. Noteworthy information is that the scaling is independent of the hardware configuration. A similar study has been realized over hyper threading. For references, the sustained speedup with 24 threads on GPP1 was 11.95 with, 13.67 for elementary contribution computation and 11.35 for the signal computation whereas those speedups are, respectively, 10.77 (+1.18x); 11.48 (+2.19x) and 10.50 (+0.85x). Hyper threading thus offers a lower benefit to the signal summation than to the rest of the implementation. Figure 14 highlights this tendency, where it is observable that the speedup is quite similar for the different configurations.

Figure 14 -Scaling on GPP1 for each configuration

Another noteworthy result from Table 3 is that the usage of MKL and ICC provide a non-negligible speed up from 1.93 to 2.52 in monothread execution. However, this speed up is down from x1.84 to x2.38 once executed over 2x6 threads on the overall execution time: some computations, like the construction of the MKL FFT plan, do not benefit as much of parallel execution.

GPU benchmarking

On the GPU side, four different GPU have been studied. Those GPU are of the same architecture, NVIDIA codename Fermi, which diverge only by their intrinsic characteristics (number of CUDA cores, core frequency, and memory bandwidth). The GPU are listed in Table 5 for reference.

In this section will be discussed the predictability of performances over a whole range of GPU for a single configuration ; then a more in depth study will focus on the impact of atomic operations over the performances. As the objective of this model is the UT field simulation resulting, the computed product is an image of such small size (100x100 float=40kb) that its transfer time is negligible before the computation time. Hence, in this section, the presented computation time do not take into account data transfer to and from the host. In Figure 15, those results are presented in graphs one for each step. It has to be noted that the normalized performances of the streaming multiprocessors of those GPUs are roughly the same for a given frequency and for a given memory bus. More precisely, this statement works best with GPU relying on the same chipset (GTX 580 and 570; GTX 480 and Tesla C2070). The redesign of the GF100 by NVidia into GF110 not only increased thermal efficiency but also benefited lightly to GPGPU performances. Moreover, this behavior is a first step toward performances prediction for a given field simulation of this algorithm on different GPU using those chipset (for those simulations), based on the measured performances for this partial set of GPUs.

Another key element to the performance is the required signal size. Experiments L0_64E and the L45 ones require a larger signal size than others to accommodate with the summation of elementary displacement over the signals (signals of size 2048 float elements versus 1024 for simpler configurations). This impacts performance at two levels: first, the size directly impacts the application of cuFFT on larger chunks of memory; secondly a wider temporal span of information means a higher collision rate when operating atomic addition to sum contributions over the signal.

To establish the impact of atomic memory operations over the simulation, Figure 16 presents the normalized execution time on a GTX 580 and a C2070 with and without them. It shows clearly that atomic addition can slow down performances, increasing execution time up to a factor 2, when memory collisions go bad. It is remarkable that performances become quasi linear to the signal size in their absence which reveals their preponderance over the rest of the computations.

Raw performances summary

In both previous sections, benchmarks have been dedicated over a single architecture. To summarize this study, Table 6 presents the best performances of this simulation over all the configurations for GPP and GPU with a temporal point of view. It appears that GPP performs better than GPU on the whole set of configurations. For modest sized problems, simulation is able to reach 10 images per second on GPP whereas GPU cap at 6.2 images per second (L0_32E). To conclude this benchmark, both implementations still require at least a speedup over an order of magnitude to reach full interactive performances. On GPP, one strategy may be to adopt a full SIMD implementation of compute-heavy sections. On GPU, work will focus on the summation of small contributions over temporal signals to reduce the influence of atomic operation.

V. Conclusion

This paper presented a new model for UT field computation, aimed at providing regularity toward massively parallel architectures. This model suffers from a lack of genericity. Indeed, the regular ray path computation requires analytical modeling of surfaces, this model is limited to canonical components, made of homogeneous and isotropic material and to direct and half-skip mode. This model has been validated, by comparison with the CIVA 11.0 software, against typical configurations and provides accurate field simulation. This new model has been implemented on both GPP and GPU, and the preliminary conclusions come as follows.

On the GPP, the implementation performances scales with the number of threads, however, this scaling is hindered on the signal formation step. By using Intel MKL and the Intel compiler, it is possible to increase the overall performances of the application by a factor 2. However, as the MKL is used for speeding up signal processing, this speedup is not shared by the other parts of the computation.

Concerning GPUs, performances of a given computation are predictable on different GPU relying on the same chipset once this configuration is benchmarked on a single hardware.

Benchmarks highlight the cost of atomic operations, which can be much slower than standard memory operations when there are a lot of collisions. Those operations amount to performances up to twice slower than standard memory operations.

Overall performances reach, for GPU, 6.2fps on modest sized configuration and can be as slow as 1.4fps on heavyweight one; whereas the GPP reaches a frame rate of 10.4fps and 3.1fps respectively. This performance is already good but still requires at least, an order of magnitude to attain interactive field simulation.

Finally, multiple improvements can be sought to speed up each parallel implementation. For example, the use of SIMD instruction on the GPP to benefit from its fine grained parallel abilities (especially on the non signal-processing steps i.e. pencil computations). Besides, by reducing the number of atomic operations in the GPU implementation, performances can improve drastically on signal summation.

Work is also in progress to provide an extensive benchmark framework, aimed to provide an auto-tuning algorithm, adapting its parameters to the class of configuration and to the hardware architecture.

Those implementations are the first step toward a fast UT simulation integrated within the CIVA 12 software.

Figure 1 -

 1 Figure 1 -Field simulation inside a weld with two different probes.

Figure 4 -

 4 Figure 4 -Formation of the impulse response at point M The impulse response is then convoluted in the frequency domain with the reference signal to obtain the resulting signal of the modulus of displacement in each point. The results of the beam formation simulation are the image of maximum of amplitude of the modulus of displacement and the image of the corresponding time of flight. Algorithm 1 presents the generic algorithm required to compute the resulting image.

Figure 5 Figure 5 -

 55 Figure5presents the algorithm of the GPP implementation, highlighting the chosen steps and the necessary call to the MKL with a specific step to predetermine the size of the resulting signals.

STARTFigure 6 -

 6 Figure 6 -Algorithm of the GPU implementationThe second steps begin with the summation of the elementary displacements to the corresponding signals for a given point P. The corresponding kernel is specialized by mode due to the differences in computations: longitudinal mode elementary responses result in real signals whereas transversal ones result in complex signals. This kernel organizes threads blocks in a similar fashion to the one evaluating the size of the resulting signals: all the threads from a block work together to compute the whole set of pencils in order to extract displacement data. Each thread then add its contribution to the corresponding displacement signals (one for each coordinate). However, as there is no way to predict the temporal span of each pencil and to avoid memory race, threads use atomic operation to contribute to the signal.

FigureFigure 8 -

 8 Figure 7 -L0-32E (left) and L45-64EP (right)

Figure 9 -

 9 Figure 9 -Temporal signal of the point with the maximum of amplitude on simulation L0-64E

  46,3/47,0 = 1,5% 10,4 mm 10,6 / 10,4 = 1,9% L45-64E 0,4 dB 20,0 mm 20,0/20,0 = 0% 2,8 mm 2,8 / 2,8 = 0,0% L45-64EP 0,4 dB 12,3 mm 12,4/12,3 = 0,8% 3,3 mm 3,3 / 3,3 = 0,0%

Figure 11 -

 11 Figure 11 -CIVA to model comparison -L45-64E In red, CIVA, in black the new model. Upper left: Simulated image -Upper right: Graph of the vertical maxima Lower left: Graph of the horizontal maxima -Lower right: Maximum signal

Figure 12 -

 12 Figure 12 -CIVA to model comparison -L45-64EP In red, CIVA, in black the new model. Upper left: Simulated image -Upper right: Graph of the vertical maxima Lower left: Graph of the horizontal maxima -Lower right: Maximum signal Figure 11 and Figure 12 highlight the focal spot measurement on the simulated image. Table2presents the measurement results for all the studied configurations and the relative error compared to CIVA. First, it is noteworthy that this model slightly overestimates the amplitude obtained by CIVA 11.0; however, in all cases, the discrepancy is less than 0.4 dB (~5%) and remains below the requirements for passing benchmarks against experimental results. Moreover, the focal spot is of the same dimension with both tools, the relative observed error is inferior to 2%. Lastly, results indicate the absence of temporal shift between CIVA and this model (as observed on previous Figures). Thus, those results validate this new UT field computation model in terms of amplitude, beam shape and time of flight.

  Figure 11 and Figure 12 highlight the focal spot measurement on the simulated image. Table2presents the measurement results for all the studied configurations and the relative error compared to CIVA. First, it is noteworthy that this model slightly overestimates the amplitude obtained by CIVA 11.0; however, in all cases, the discrepancy is less than 0.4 dB (~5%) and remains below the requirements for passing benchmarks against experimental results. Moreover, the focal spot is of the same dimension with both tools, the relative observed error is inferior to 2%. Lastly, results indicate the absence of temporal shift between CIVA and this model (as observed on previous Figures). Thus, those results validate this new UT field computation model in terms of amplitude, beam shape and time of flight.

Figure 15 -

 15 Figure 15 -Graphical representations of GPU normalized execution time This graph presents signal computation time normalized relatively to the signal size

Figure 16 -

 16 Figure 16 -GPU normalized time with and without atomic operations (C2070 and GTX580 only) This graph presents signal computation time normalized relatively to the signal size

Table 3 -GPP Results -speedup for ICC & MKL -hardware GPP1 -2x Xeon 5590 6 cores without HT

 3

Table 4

 4

	presents the time repartition over executing GPP
	implementations on L0_64E and L45_64E. This confirms
	that the MKL execution benefit greatly to the impulse
	response signal manipulations; however the other
	computations (generating elementarily impulse response,
	and signal size determination) do not benefits as much from
	the MKL as their overall proportions increase.

Table 4 -Execution time repartition on GPP1

 4

Table 5 -GPUs studied

 5

		Chip	# of SM	# of cores	Core frequency	Memory bandwidth
	GTX 480	GF100	15	480 1215Mhz 177.4GB/s
	Tesla C2070	GF100	14	448 1150Mhz 144GB/s
	GTX 570	GF110	15	480 1464Mhz 152 GB/s
	GTX 580	GF110	16	512 1554Mhz	192.4 GB/s

Table 6 -Performances summary (GPU and GPP)

 6

	GPU using GTX 580;
	GPP using GPP1 -2x GPP Intel Xeon 5590 with Intel ICC and MKL;

MSVC Version 16.00.40219.01 with Visual Studio

SP1 2 Intel(R) C++ Intel(R) 64 Compiler XE, Version 14.0.0.103 Build 20130728