C Kervazo

T Liaudat

Jerome Bobin

Faster and better sparse blind source separation through mini-batch optimization

Keywords: Blind source separation, Sparse representations, Alternating least-squares, Mini-batches, Matrix factorization, Robust estimator aggregation, Riemannian manifold

Faster and better sparse blind source separation through mini-batch optimization

C Kervazo, T Liaudat, Jerome Bobin

Introduction

Towards large-scale BSS

During the last decades, Blind source separation (BSS) has become a major analysis tool to learn meaningful decompositions of multivalued data, in a wide variety of scientific fields such as audio processing [START_REF] Vincent | Probabilistic modeling paradigms for audio source separation[END_REF][START_REF] Févotte | Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis[END_REF], biomedical data processing [START_REF] Negro | Multichannel intramuscular and surface EMG decomposition by convolutive blind source separation[END_REF] or astrophysics [START_REF] Bobin | Joint planck and wmap cmb map reconstruction[END_REF], to only cite three of them. According to the standard linear mixture model, the data are composed of m observations, which stem from linear combinations of n elementary signals called sources, each of them having t samples. In matrix form, the model then writes as:

X = AS + N (1)
where X (size m × t) is called the observation matrix and N is the noise matrix accounting for model perturbations. As such, the goal of BSS is to retrieve the mixing coefficients A (size m × n) as well as the sources S (size n × t) up to a mere scaling and permutation indeterminacy. BSS is a special case of constrained matrix factorization. Being an ill-posed problem, extra prior information needs to be exploited, such as the statistical independence of the sources (Independent Component Analysis family -ICA [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]), the nonnegativity of A and S (Non-negative Matrix Factorization -NMF [START_REF] Gillis | Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization[END_REF][START_REF] Gillis | Successive nonnegative projection algorithm for robust nonnegative blind source separation[END_REF][START_REF] Gillis | The why and how of nonnegative matrix factorization[END_REF][START_REF] Vandaele | Heuristics for exact nonnegative matrix factorization[END_REF]).

In this work, we will focus on sparse modelling [START_REF] Zibulevsky | Blind source separation by sparse decomposition in a signal dictionary[END_REF][START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF][START_REF] Le Roux | Sparse nmf-half-baked or well done?[END_REF] and assume the sources to be sparse in some domain Φ S . This signal representation should be adapted to the geometrical content of the sources to be separated [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF].

The rapid increase of the data size to be analysed mandates the development of fast and scalable algorithms. So far, sparse BSS has mainly been applied to small (e.g the space telescope Chandra [START_REF] Picquenot | Novel method for component separation of extended sources in x-ray astronomy[END_REF] -a few million pixels) to middle size data (e.g the data of the space mission Planck [START_REF] Bobin | Joint planck and wmap cmb map reconstruction[END_REF] -half a billion pixels).

As such, although currently available methods can be very effective, they are not suited to cope with the data challenges to come in this field (e.g the Square Kilometer Array radio-telescope [START_REF] Chapman | The scale of the problem: recovering images of reionization with generalized morphological component analysis[END_REF] will provide several billions of pixels). The main challenge is then how to design a fast and distributed sparse BSS algorithm without sacrificing the separation performances ?

BSS as an unsupervised matrix factorization problem

In this article, we will focus on tackling noisy sparse BSS problems. For that purpose, the following standard assumptions are made:

• The noise is assumed be Gaussian additive independently and identically distributed 1 . This naturally entails that the cost function will first be composed of a quadratic data fidelity term: 1 2 X -AS 2 F , where the Frobenius norm . F is used to measure the discrepancy between the data and the model.

• The sources are assumed to be sparse in some signal representation Φ S (matrix of size T × t, with T t). For the sake of simplicity, Φ S will be assumed to be an orthogonal representation. This could also be relaxed to more general tight frames [START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF]. The sparsity of the sources is classically measured by a re-weighted 1 -norm : R S (SΦ T S ) 1 . The operator the denotes the Hadamard product. The regularization parameters R S (size n × T ) control the trade-off between the data fidelity and the sparsity terms.

• To alleviate degenerated solutions ( A F → ∞ and S F → 0) due the usual scale indeterminacy in BSS, the last term enforces the mixing matrix to belong to the oblique ensemble: all the columns A i , i ∈ [1, n] of the mixing matrix must lie on the 2 unit hypersphere. The characteristic function of an arbitrary set U is denoted as ι U (.).

Altogether, these assumptions lead to the following generic formulation of the sparse BSS problem:

minimize A∈R m×n ,S∈R n×t 1 2 X -AS 2 F + R S (SΦ T S ) 1 + ι {∀i∈[1,n], A i 2 2 =1} (A), (2)
Being a multi-convex matrix factorization [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF] problem, it is generally tackled by sequentially alternating minimization steps with respect to the mixing matrix A and the sources S.

Beyond sparse BSS, a standard strategy to cope with large-scale data in generic matrix factorization problems of the form (2) consists in processing small size batches (i.e. mini-batches) in a distributed manner. As such, the data are decomposed into sub-matrices (sets of columns of X). So far, the generic distributed sparse matrix factorization algorithms in the literature can be decomposed into two categories :

• Stochastic approximations, which have been introduced in the general framework of dictionary learning. In [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF], the authors proposed an online algorithm in which the dictionary (i.e. equivalently the mixing matrix) is computed only by minimizing an upperbound of the empirical cost. While the online setting is a special case with mini-batches size of t b = 1, the algorithm is generalized to arbitrary values of t b . More recently, this work has been extended in [START_REF] Mensch | Stochastic subsampling for factorizing huge matrices[END_REF] to cope with datasets that, in our context, would be large both in the number of samples t and number of observations m.

• Stochastic gradient descent -SGD [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF]: [START_REF] Davis | The sound of apalm clapping: Faster nonsmooth nonconvex optimization with stochastic asynchronous palm[END_REF] extended the use of the Proximal Alternating Linearized Minimization (PALM) to mini-batches, making it possible to tackle large datasets. Such an approach has also been extended to the case of hyperspectral imaging in [START_REF] Thouvenin | Partially asynchronous distributed unmixing of hyperspectral images[END_REF], in which the authors argue a higher flexibility than in [START_REF] Davis | The sound of apalm clapping: Faster nonsmooth nonconvex optimization with stochastic asynchronous palm[END_REF] (where the asynchronicity has a high impact on the allowable step sizes, counter-balancing its positive effects).

On the one hand, these generic algorithms can virtually tackle any sparsityregularized matrix factorization problem with a quadratic data fidelity term.

They have however not been particularly adapted for sBSS problems, where the optimization strategy plays a key role to provide satisfactory results [START_REF] Kervazo | Use of palm for 1 sparse matrix factorization: Difficulty and rationalization of an heuristic approach[END_REF].

Furthermore, they are generally based on gradient descent, for which such a tuning is cumbersome [START_REF] Kervazo | Use of palm for 1 sparse matrix factorization: Difficulty and rationalization of an heuristic approach[END_REF].

On the other hand, the GMCA algorithm [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF] is a dedicated sparse BSS algorithm. In contrast, it is based on a projected Alternating Least Square (pALS) scheme, which originates from early works in NMF [START_REF] Paatero | Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values[END_REF]. More precisely, when updating one of the variables A or S, a full minimization of the (differentiable) data-fidelity term is carried out, and then a projection over the non-differentiable constraints is performed by applying the corresponding proximal operator [START_REF] Parikh | Proximal algorithms[END_REF]. In this context, the GMCA has been quite successful as it comes with an automatic scheme to fix the regularization parameters R S . This procedure has been showed to increase its robustness to initialization, local minima and noise. Unfortunately, the GMCA algorithm does not computationally scale very well to large size data.

As such, existing state of the art methods do not currently handle the largescale sBSS problem properly: on the one hand classical sBSS methods do not scale well with the data size, while generic matrix factorization methods, although computationally efficient, do not usually gives good results in the context of sBSS.

Contributions

The main contribution of this article is to introduce a new algorithm coined distributed GMCA (dGMCA), which combines projected alternate least-squares with mini-batch optimization. The dGMCA algorithm builds upon the parallel computation of estimates of the mixing matrix A from data mini-batches. The originality of the dGMCA algorithm then lies in the construction of an aggregated estimate using a robust mean on the hypersphere, which better respects the Riemannian geometry of the oblique constraint. Numerical experiments have been carried out on various mixing scenarios. They first highlight that the proposed algorithm provides a scalable sparse BSS algorithm without sacrificing the quality of the factorization.

Furthermore, and more surprisingly, the dGMCA algorithm can significantly outperform standard sparse BSS algorithm in challenging problems, such as cases where the sources have a highly sparse representation. We discuss the connections between dGMCA and SGD, which gives some insight into the behavior of dGMCA based on the exploratory power of stochastic minibatch optimization and the resulting implicit regularization. The dGMCA algorithm is further illustrated on realistic γ-spectroscopy data. Preliminary results were presented in the SPARS 2019 conference [START_REF] Liaudat | Distributed sparse bss for large-scale datasets[END_REF].

Notations

In the following, scalars will be written in lower case letters a, vectors in bold a and matrices as A. A j stands for the j th column of A, while A i is the i th line. These notation are extended to subsets of lines or columns:

if J b is a subset of [1, t], S J b denotes the columns indexed by J b . Generally speaking, we will consider B subsets J b , b ∈ [1, B] that form a partition of [1, t], which will be used to index mini-batches. We will furthermore write #J b the number of elements in J b , that is in each mini-batch.

For each mini-batch, we will have access to several estimations of the same matrix A. Each one will be denoted by

A[J b ], b ∈ [1, B]. The aggregated estimate of the different A J b is denoted as Â.
In iterative algorithm, we will write a (k) the estimation at iteration k of the variable a.

2. dGMCA: a distributed sparse alternating least-squares algorithm

The distributed GMCA algorithm

When it comes to distributed computation, a straightforward approach consists in performing the most intensive computational burden on smaller chunks of data. Following generic distributed matrix factorization algorithms, the proposed method performs by first splitting the data matrix X into B disjoint sub-matrices. Without loss of generality, the number of samples t is assumed to be a multiple of B, thus leading to a complete decomposition of the full data X with constant batch size. Therefore, each of the mini-batches has

t b = t/B columns (#J b = t b ).
In the GMCA algorithm, the updates of mixing matrix A and the source matrix S have the largest computational cost. Each iteration k can be described by the following two steps:

1 -S is updated assuming a fixed A.

S (k) = S R S A (k-1) † X , (3)
where A (k-1) † is the pseudo-inverse of A (k-1) . The operator S R S is the soft-thresholding operator with thresholding parameters R S .

2 -A is updated assuming a fixed S:

A (k) = Π . 2 =1 XS (k) † , (4)
where Π . 2 =1 is the projection onto the m-sphere.

In the proposed dGMCA, each main iteration or epoch 2 of the distributed GMCA algorithm can be decomposed into the following two stages:

• The first stage amounts to performing an estimation of the sources and the mixing matrix independently from each batch

X J b , b ∈ [1, B]
. This 2 If one uses the standard machine learning vocabulary.

is the most computationally intensive part of the algorithm and it is performed in a distributed way.

• Each mini-batch leads to an independent estimate {A[J b ]} b∈ [1,B] of the mixing matrix. The second stage consists in combining or aggregating these different estimates to eventually produce a single one.

More formally, the dGMCA algorithm is summarized in Algorithm 1 The Algorithm 1 dGMCA

1: for k = 1, • • • , K do 2: Choose J 1 , J 2 , • • • , J B as a partition of [1, t] 3: for b = 1, • • • , B do 4: ŜJ b (k) = S R S (k) ( Â(k-1) † X J b (k) )
Use Eq. ( 12) for R S choice

5: A[J b (k)] = Π . 2 =1 (X J b (k) ŜJ b (k) † ) 6:
end for 7:

Â(k) = AGGREGATE(A[J 1 (k)], A[J 2 (k)], ..., A[J B (k)])
8: end for

9: return Â(K) , Ŝ(K)
function AGGREGATE is used as a generic terminology for the aggregation procedure, which will be described in details in the following section.

Manifold-based mixing matrix aggregation

The objective of the aggregation step is to build a single estimate of A from the various estimates

A[J b ], b ∈ [1, B].
For that purpose, the aggregated estimate  can naturally be defined as the barycenter of the different estimates according to some distance φ, with some weights {ω b } b∈ [1,b] as follows:

 = Argmin A B b=1 ω b φ (A, A[J b ]) , (5)
where the barycentric weights are positive and sum to one

: ∀b ∈ [1, b], ω b ≥ 0 and B b=1 ω b = 1. A straightforward choice for φ could be the standard Euclidean distance : φ (A, A[J b ]) = n j=1 A j -A j [J b ] 2
2 , which eventually defines  as a standard weighted average of the estimates. However, a key property of ( 2) is that the mixing matrix should belong to the Oblique manifold, which would not be necessarily satisfied with a standard Euclidean barycenter. For that purpose, dedicated distances will be rather used to precisely take into account the Riemannian geometry of the Oblique constraint.

Frechet mean on the hypersphere

The oblique constraint implies that each column of the mixing matrix A belongs to the hypersphere of dimension m -1 or m -1-sphere S m-1 , which is a Riemannian manifold (see [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] for more details). Building an aggregation procedure that respects the underlying Riemannian geometry is naturally done by choosing φ as the geodesic distance on the m -1-sphere.

The extension of the notion of barycenter of data points onto a Riemannian manifold equipped with its geodesic distance is known as the Fréchet mean.

Precisely, for each column j of the mixing matrix, the Fréchet mean of the

individual estimators A[J b ] is defined as: ∀j ∈ [1, n], Âj = Argmin a∈R m b∈[1,B] ω b d β a, A[J b ] j , (6)
where 1 ≤ β < +∞. The case β = 2 corresponds to the 2 norm along geodesics of the m -1-sphere.

Finding the solution to the above minimization problem is not straightforward but can be computed thanks to a dedicated iterative gradient descent algorithm [START_REF] Asfari | On the convergence of gradient descent for finding the riemannian center of mass[END_REF]. Under some conditions, this algorithm is proved to converge to a critical point of the above problem. More precisely, each iteration of Afari's algorithm (cf. Algorithm 2) are decomposed into two stages:

• Computation of the gradient :

for β ≥ 1 the gradient ∇J β ( Âj (k) ) of the cost function J β ( Âj (k) ) = B b=1 ω b d β Âj (k) , A j [J b ] is defined as : ∇J β ( Âj (k) ) = - B b=1 ω b d β-2 Âj (k) , A j [J b ] log Âj (k) A j [J b ] , (7)
where the logmap log Âj (k) is roughly speaking the projection onto the tangent plane of the m -1-sphere about Âj (k) [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF].

• Update of the current mean estimate : the updated estimate is then defined as the exponential map about A j [J b ] applied to the gradient with step size ρ3 . This step merely "back-projects" the gradient -which belongs to the tangent space about Âj (k) -onto the m -1-sphere.

Robust Fréchet mean on the hypersphere

The use of small mini-batches makes the separation process more prone to generate outliers in the estimated A[J b ] (cf. Section 3.4 and numerical experiments of Section 3). Unfortunately, the Fréchet mean is not robust to such outliers. As such, more robust distances have been used [START_REF] Fletcher | Robust statistics on riemannian manifolds via the geometric median[END_REF][START_REF] Arnaudon | Medians and means in riemannian geometry: existence, uniqueness and computation[END_REF].

In the present setting, a natural choice would be to choose β = 1, which

for j = 1, • • • , n do 3:
while convergence is not reached do 4:

∇J β ( Âj (k) ) = -B b=1 ω b log Âj (k) (A j [J b ])
5:

Âj (k+1) = exp Âj (k) -ρ∇J β ( Âj (k) )) 6:
end while

7:
end for 8: end procedure corresponds to the usual 1 norm. However, the 1 norm is not differentiable about 0. In dGMCA, this is mitigated by building a smooth approximation d 1 µ of d 1 based on Nesterov's smoothing technique [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]:

d 1 µ (a, b) = Argmax u ∞≤1 < a -b, u > - µ 2 u 2 2 . ( 8
)
For a, b ∈ R m . This approximated distance is differentiable and its gradient is µ-Lipschitz4 . This entails that the gradient of the cost function takes the form:

∇J β ( Âj (k) ) = - B b=1 ω p ∇d 1 µ Âj (k) -A j [J b ] log Âj (k) A j [J b ] , (9)
where the gradient of d 1 µ is given by [START_REF] Becker | Nesta: a fast and accurate first-order method for sparse recovery[END_REF] for all i ∈ [1, m]:

∇d 1 µ ( Âj (k) ) i =      µ -1 Âj (k) i , if | Âj (k) i | < µ, sign( Âj (k) i ), otherwise. (10)
The resulting algorithm is similar to Algorithm 2, with the exception that the distance and its gradient are evaluated according to Equations 8 and 10.

Implementation details

Use of a reweighted aggregation

Using mini-batches of small size generally yields very diverse estimates of the mixing matrix. For instance, in mini-batches where a source is dimly active or inactive at all, the estimated mixing matrix is likely to exhibit a degraded condition number, which eventually leads to a dramatic increase of the noise that contaminates the estimated sources. Therefore, the barycentric weights used for either the Fréchet mean or its robust version are proportional to the noise variance of the estimated sources. Let's denote by

W[J b ] = A[J b ]
+ the pseudo-inverse of the mixing matrix estimated from batch J b , each weight is defined as ω j b :

∀b ∈ [1, B]; ω j b = W[J b ] j Σ N W[J b ] j T -1 B p=1 W[J p ] j Σ N W[J p ] j T -1 (11)
where Σ N is the data noise covariance matrix and

W[J b ] j is the j-th line of W[J b ].
The maximum number of iterations of the (robust) Fréchet mean evaluation algorithm was fixed to 1000 with a stopping rule so that the algorithm stops when the angle between two iterates is lower than 10 -6 .

Choosing the regularization parameters

As pointed out earlier, the regularization parameters play a key role, as they correspond to thresholds applied to the sources. Thus, they select specific entries of the sources during the separation process. Following the morphological diversity principle [START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF], the last versions of the GMCA algorithm make use of a strategy where the percentile of retained entries of S (k) linearly increases from step to step, which enables to keep an increasing number of entries in a way which is well adapted to the actual distribution of the sources. This strategy has been shown to lead to a high robustness of the algorithm with respect to spurious local critical points. However, to be applied it requires the knowledge of all the source samples, which is impossible in the proposed distributed optimization procedure, where the mini-batches are processed separately (i.e. without communication between the distributed processes). To that end, a strategy based on an exponential decay of the parameters is used:

R i(k) i S i = κσ i + Ŝ(k) i ∞ -κσ i exp(-kα i ), (12)
where σ i is the estimated standard deviation5 of the noise contaminating the source S i , the constant κ is generally chosen between 1 and 3 based on a desired probability of false discoveries and α i is a parameter controlling the exponential decay decrease. In practice, α i is fixed to 2 for all the sources; other choices close to the unit do not change significantly the separation results. The proposed strategy is adapted to distributed computing as

Ŝ(k) i ∞
can be evaluated as the maximum over each mini-batch:

Ŝ(k) i ∞ = max b∈[1,B] ŜJ b (k) i ∞
. As well, σ i can be evaluated empirically as the median of its estimations over the mini-batches:

σ i = median b∈[1,b] σ i [J b (k)].
dGMCA convergence and stopping rule

From the theoretical point of view, the convergence of the distributed GMCA can hardly be proved since: i) it is based on an inexact minimization scheme that alternates least-squares estimates and projections/proximal step, and ii) for the sake of robustness, the thresholding strategy is based on the decrease of the regularization parameter value during the separation process, and iii) each iterate is obtained as a non-linear aggregation of batch-based estimates. However, experiments tend to show that the algorithm stabilizes after a certain number of iterates. If convergence cannot be rigorously claimed, there are empirical clues that support the stabilization of the algorithm in practice. To obtain convergence guarantees, one could envision to incorporate dGMCA as a warm-up step of a two-step algorithm, similarly as in [START_REF] Kervazo | Use of palm for 1 sparse matrix factorization: Difficulty and rationalization of an heuristic approach[END_REF]. The second refinement step would then be an asynchronous PALM algorithm [START_REF] Davis | The sound of apalm clapping: Faster nonsmooth nonconvex optimization with stochastic asynchronous palm[END_REF]. This is out of the scope of this article and left to future work.

The stopping rule is based on the maximum (across columns) angular distance between two iterates of the estimated mixing matrix. The algorithm stops when such distance is smaller than 10 -6 . The maximum number of iterations is fixed to 10 4 . As a pre-processing, the columns of the data may be randomized to provide more homogeneous batches. Similarly, randomization has also been tested at each iteration to mimic stochastic mini-batch optimization in the projected ALS framework. For the experiments below, randomisation at each iteration did not provide further significant improvements.

Numerical experiments

In this section, the performances of the dGMCA algorithm are investigated in various experimental settings, on both simulated and realistic data.

Experiments on simulated data

Comparison set-up

In this subsection, we first make use of synthetic random data, which allows to assess the robustness of the different methods in various experimental scenarios. The influence of several parameters has been evaluated:

the number of sources n and observations m, the condition number of the mixing matrix and the sparsity level of the sources. Since all led to similar conclusions, we hereafter focus more specifically on the two last items, which yield to the most insightful results. To that end, the data are synthesized as follows:

• The entries of the sources S j , j ∈ [1, n] are independently and identically distributed according to a Generalized Gaussian distribution with parameter 0 < γ ≤ 1.

• The mixing matrix is picked at random from a Gaussian distribution, and further processed to have columns with unit 2 norm and a predefined condition number.

Unless stated differently, each single experimental result will be given as the mean over 25 Monte-Carlo simulations with different mixing matrices, sources and noise realizations.

Comparisons will be carried out with the following unsupervised matrix factorization algorithms:

• GMCA : This non-distributed algorithm [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF] serves as a baseline to evaluate the proposed distributed separation procedure.

• Online dictionary learning : This algorithm [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF] is a classical one for solving large-scale sparse matrix factorization problems. For fairer comparisons, the regularization parameter has been optimised based on several simulations.

• distributed GMCA : the dGMCA will come with two distinct aggregation procedures, namely with the Fréchet mean and its robust alternative.

To assess the separation quality, the mixing matrix criterion C A is used:

C A = mean(|P Â † A -I|) (13)
With A the true mixing matrix and P Â † the pseudo-inverse of the solution found by the algorithm corrected through P for the scale and permutation indeterminacy. The mean is the average of all the elements inside the matrix.

In contrast to more standard criteria that are based on the estimated sources, the advantage of the mixing matrix criterion is that it is less sensitive to the regularization of the sources, which generally differs between algorithms.

Sparsity level ρ

The sparsity level ρ of the sources impacts the separation process in two ways:

• The sources are generated so as to be statistically stationary. However, in the very sparse regime (e.g. ρ = 0.1), the values taken by a single small mini-batch may largely change between two realizations. This can lead to outlier mini-batches that will impact the aggregation process. In contrast, mildly sparse sources should lead to a more stable aggregation procedure;

• Very sparse sources tend to lead to sharper critical points that are much harder to escape. These local minima are more likely to be smoothed out when the sparsity level decreases.

The number of observations is set to 20, the number of sources to 5 and the number of samples per sources to 10000. The signal-to-noise ratio is fixed to 40dB. The condition number of the mixing matrix is set to 3. As an illustration, Figure 1 displays the evolution of the mixing matrix criterion These results reveal two first distinct regimes: i) For mildly sparse sources (left panel of Fig. 1), the GMCA algorithm yields better results and there is only a slight discrepancy between the two dGMCA methods. This highlights that equipped with the (robust) Fréchet mean, combining mini-batch optimization and aggregation leads no performance loss: dGMCA allows to perform distributed computation without deteriorating significantly the separation quality.

ii) For less sparse sources (cf. right panel of Fig. 1), the results of robust dGMCA are consistent with the previous regime: robust dGMCA reaches a separation quality that is very close to GMCA for even smaller mini-batches. One of the differences with regime i) is that GMCA obtains slightly worse results. This is expected, as more partial correlations may occur in this regime. Furthermore, the discrepancy between Fréchet mean aggregation and its robust counterpart increases again strongly when the batch size becomes smaller than 100 samples. This might also be due to an increased number of partial correlations: small mini-batches containing such samples might be much more difficult to unmix, generating more inhomogeneous mini-batches.

To go further, Fig. 2 shows the evolution of the mixing matrix criterion as a function of the sparsity level for two different mini-batch sizes: t b = 10 and t b = 100 samples. For ρ > 0.2, this figure confirms the above comments:

equipped with the robust Fréchet mean, the dGMCA and GMCA algorithms provide very close results. As well, the proposed robust aggregation procedure yields significantly better separation quality.

When the sources are highly sparse (typically for ρ < 0.1), the performances of the GMCA algorithm rapidly degrade. Very likely, this originates from the sharpness of the critical points, which are much harder to escape. Quite astonishingly, in this precise case, the dGMCA algorithm with the robust Fréchet mean aggregation performs much better than GMCA for small batch size. This surprising results originate from the exploratory power of small batch sizes, which better prevents the dGMCA algorithm to be trapped in sharp critical points. This phenomenon will be discussed in more details in Section 3.4.

Condition number

Ill-conditioned mixtures, as measured by the condition number of the mixing matrix A, generally lead to arduous separation problems. Mixing matrices with large condition numbers yield two major bottlenecks: i) an increased noise level in the source domain, and ii) the mixtures are closer to co-linearity, which makes the sources harder to distinguish.

In these experiments, the noise level is fixed to 40 dB and the sparsity level is ρ = 0.1. The number of observations is set to 20, the number of sources to 5 and the number of samples per sources to 10000. Figure 3 shows the evolution of the mixing matrix criterion C A as a function of the mini-batch size t b for two values of the condition numbers: left panel 2.5 and right panel 7. As expected, the quality of the separation results of all methods decrease when the condition number increases. Similarly to the tests performed in the previous section, the dGMCA algorithm has better results for relatively small mini-batch sizes (but when the Fréchet mean is used, it eventually deteriorates for t b < 25, cf. Fig. 3). The use of small batches along with the robust Fréchet mean leads to an improvement for t b < 25, which becomes more significant when the condition number increases up to a gain of about one order of magnitude. Similarly, when the mini-batch size decreases, the discrepancy between the two methods increases.

Implementation and computation time

We now conclude this experimental section with the computation time of dGMCA, which depends both on the complexity of one epoch and the number of required epochs.

O(b(mn + n 2 + m) + t b nmK), (14)
where the last term correspond to the Fréchet mean and K corresponds to the number of iterations required for its computation. As such, except for very small mini-batches, the linear gain of using dGMCA over GMCA dominates. This is experimentally confirmed (cf. Fig. 4): in practice the computation time for a given number of iteration does not deviate much from linearity (in particular, the transfer costs between the nodes are negligible in comparison to the computation time).

Application to γ-ray spectroscopy realistic simulations: sparse case

In this section, the behavior of the dGMCA algorithm is evaluated in the context of γ-ray spectroscopy. This is one of the main methods used for measuring the activity concentrations of radionuclides in environmental samples. It particularly plays a key role to monitor the radiological environment or perform radioecology studies and nuclear incident preparedness. A γ-ray spectrum is the histogram of the number of detected γ-ray photons in the detectors. In this context, an observation is formed by the linear combination of the contribution from various radionuclides. Each one is described by a signature in energy that is composed of one or several emission lines to which a compton continuum is associated as displayed in Figure 5. The goal of this experiment is to jointly estimate the activity of each radionuclide (a.k.a. the mixing matrix) as well as their signature (a.k.a. the sources) from several observations. These simulations are composed of 5 radionuclides: 7 Be, 22 Na, 40 K, 137 Cs, 210 Pb, which are representative of aerosol samples [START_REF] Xu | Sparse spectral unmixing for activity estimation in γ-ray spectrometry applied to environmental measurements[END_REF] and featured in Figure 5. The number of observations is fixed to m = 20 and the number of samples per source is equal to 16940.

These data are particularly interesting as they allow to evaluate the performances of dGMCA when the samples are clearly non-stationary, highly sparse and with a large dynamic range (the source samples basically span 2 to 3 orders of magnitude). The sources are modeled in the wavelet domain:

γ-ray observations are first decomposed into an undecimated uni-dimensional wavelet frame [START_REF] Starck | The undecimated wavelet decomposition and its reconstruction[END_REF] before applying any BSS method. The number of scales is fixed to 5, which yields a number of wavelet coefficients that is equal to 81200; these are obviously not large-scale data but it already allows to highlight some remarkable results.

Figure 6 shows the reconstructed solution with GMCA, ODL and dGMCA equipped with the robust Fréchet mean with t b = 10; it also displays the estimation error in transparent solide line. This figure first shows that dGMCA provides a very good reconstruction of 2 2 Na signature, while both GMCA and ODL exhibit clear leakage from other sources. The estimation error of the dGMCA solution does not present any structure and is mainly dominated by noise.

Figure 7 features the evolution of the mixing matrix criterion as a function of the mini-batch size for two different levels of the signal-to-noise ratio: 40 and 80 dB. These values might seem large but it has to be recalled that the dynamic range is very large; a small amount of noise might already erase a significant part of the compton continuum while leaving only the photon peaks. This experiment first shows that GMCA, ODL and dGMCA with the standard Fréchet mean performs rather poorly. In agreement with the results of the previous subsection studying the impact of the sparsity level ρ, the use of the robust aggregation makes the dGMCA algorithm largely outperforms these methods, especially when the mini-batch size is smaller than t b = 50. Further randomizing the mini-batches entails an extra improvement, especially for middle-size mini-batches for 100 < t b < 1000. The gain is particularly large when the noise level is small.

To explain the results of such a setting, we advocate the fact that the sources are by a large extent dominated by few photon peaks, which is likely to create spurious critical points that are hard to escape. This might explain why neither the GMCA algorithm nor the dGMCA algorithm without the robust aggregation are able to perform correctly. The use of the robust Fréchet mean makes the algorithm less prone to be stuck in sharp critical points.

This will be discussed in the next section. From these experiments, two distinct regimes can be highlighted:

• Mildly sparse sources: in this setting, the robust dGMCA and GMCA algorithms perform similarly; going distributed comes at almost no cost as soon as the mini-batch size is large enough. For small mini-batch sizes (typically smaller than t b = 100 and even for very small mini-batch sizes), the dGMCA algorithm with the robust aggregation leads to performances that are very close to GMCA. This shows that the proposed dGMCA is an efficient approach for separation of large-scale mildly sparse signals, which are typical of natural images with a multiscale representation such as wavelets.

• Sparse to very sparse sources: As expected, small mini-batch sizes lead to more heterogeneous mini-batches, which are more likely to be seen as outliers. In this case, the proposed robust Fréchet mean provides superior and more robust separation performances with respect to the standard Fréchet mean. Interestingly, using small mini-batch sizes with dGMCA in its robust version leads to significantly enhanced results with respect to the GMCA algorithm, which can also be observed for ill-conditioned mixtures. We now give insights concerning this astonishing phenomenon.

Stochasticity of mini-batch optimization and implicit regularization

Connections with mini-batch GD :. A first remark is that either for the update of the mixing matrix or for the sources -and ignoring the projection step -a least-square update is virtually equivalent to a Newton iteration with gradient path length equal to 1. As such, our explanation will be based on links between Gradient Descent (GD) and ALS. More precisely, mini-batch optimization has been a widespread strategy in minimization methods relying on GD6 .

Mini-batch GD induces implicit regularization :. Understanding the impact of optimization in learning parameters from matrix factorization problems [START_REF] Gunasekar | Implicit regularization in matrix factorization[END_REF] or (deep) neural networks [START_REF] Neyshabur | Implicit regularization in deep learning[END_REF] has attracted a lot of interest during the last few years. More specifically, it has been emphasized that a specific optimization method might enable some implicit regularization, which tends to favor the convergence towards critical points with specific properties. Such a phenomenon is well known in the context of machine learning, where it has been noticed for long that using small stochastic mini-batches along with GD can improve the results over full batch methods [START_REF] Lecun | Efficient backprop[END_REF][START_REF] Keskar | On large-batch training for deep learning: Generalization gap and sharp minima[END_REF][START_REF] Hardt | Train faster, generalize better: Stability of stochastic gradient descent[END_REF]. A common interpretation is that using small size mini-batches is important as it injects noise in the optimization process, which is essential to escape certain types of critical points. More specifically, it has been empirically shown in [START_REF] Xing | A walk with sgd[END_REF] that using stochastic mini-batches enables to explore broader areas further away from the initial point. The authors also argue that the stochasticity introduced by small size mini-batches through a structured noise strongly favors flat minimizers that are akin to generalize better [START_REF] Hochreiter | Flat minima[END_REF].

As for BSS, the factorization problem at play is furthermore determined or over-determined, which means that the optimization landscape is likely to be largely different and probably less smooth, leading to the presence of spurious local critical points. In this setting, regularization, either being explicit or implicit, is of paramount importance. To that respect, the case of highly sparse sources is particularly illustrative (e.g. the γ-spectroscopy application): dGMCA performs very well as the mini-batch size decreases, while GMCA does not. As such, a similar implicit regularization is very likely to be at play within dGMCA, for which the generalization notion would translate into minimizers that are less sensitive to a given realization of the sources. Indeed, as shown by the γ-ray spectroscopy example, the algorithm performs much better in the highly sparse case, where spurious critical points tend to be sharper and more difficult to escape. The stochasticity induced by mini-batches is likely to prevent the dGMCA algorithm to be stuck in such solutions of the optimization landscape.

To further illustrate this phenomenon, Figure 8 displays the histograms of the mixing matrix criterion over the different mini-batches after 1000 iterations of dGMCA (close to "convergence") when the Fréchet mean (left panel) or its robust version (right panel) is used. In this experiment, the number of observations is set to 20, the number of sources to 5, the sparsity level to ρ = 0.1, and the condition number of mixing matrix to 1. The number of samples per source is fixed to In the first case, one can notice that using large size mini-batches provide more stable solutions, with a smaller scatter of the criterion across mini-batches. In contrast, using smaller minibatches leads to a broader exploration of the parameter space as testified by more widespread values of the mixing matrix criterion. As most of the values are quite poor (close to 10 dB), the aggregated estimate is not satisfactory.

Switching to robust aggregation (right panel Figure 8), the observed scatter is very similar but now the aggregated value is much more robust to outlier mini-batches. At first sight, it might look strange to produce an aggregated estimate that is much better than the majority of each individual estimate: it is largely on the left side of the distribution. However, it is likely that the exploration of the optimization landscape goes in random directions, which are not measured in this histogram. It further reveals that the aggregation step is essential to capture an average (and robust) estimate out of this stochastic exploration.

Figure 9 shows the evolution of the mixing matrix criterion when the minibatch size increases from 10 to 50. 10 3 for B = 10, 100 and 1000. One can first point out that for large mini-batch size (i.e. t b = 1000), the separation accuracy is poor and does not improve when the number of mini-batches increases. In this regime, it is likely that the various estimates of the mixing matrix do not present enough stochasticity to prevent the algorithm from being stuck in a spurious critical point. For a larger number of mini-batches of smaller sizes, the exploration power highly increases. This define a clearly distinct regime, which suggests some "phase transition" with respect to the mini-batch size for a fixed number of mini-batches.

To sum up, an optimization landscape exploratory phenomenon related to the stochasticity of mini-batch optimization occurs when small mini-batches are considered (typically a few times the number of sources n). In this regime, the separation quality will improve when the number of mini-batches increases. This phenomenon will vanish when the mini-batch size is too large.

Software

The dGMCA algorithm will be made available online as part of the pyGM-CALab toolbox (https://github.com/jbobin/pyGMCALab).

C

  A with respect to the mini-batch size for two values of the sparsity level of the mixing matrix: ρ = 0.25 in the left panel and ρ = 0.5 in the right panel.

5 Figure 1 :Figure 2 :

 512 Figure 1: Evolution of the mixing matrix criterion as a function of the sparsity level with ρ = 0.25 (left panel) and ρ = 0.5.

Figure 3 :

 3 Figure 3: Evolution of the mixing matrix criterion as a function of the mini-batch size for two distinct values of the mixing matrix condition number.

Figure 4 :

 4 Figure 4: Computation time gain of the dGMCA algorithm with respect to the GMCA algorithm as a function of the reduction gain t/B. The dGMCA algorithm has been run on a PC equipped with 8 Amd CPUs, each one has 6 cores Istanbul Opteron 8431 at 2, 4Ghz. The parallelization of the code has been carried out via python/C++ wrappers with OpenMP.

3. 4 .

 4 Discussion -robustness and implicit regularization 3.4.1. Summing-up the experimental results

Figure 5 :

 5 Figure 5: γ-ray spectroscopy: example of a single observation and the contribution of each of the radionuclide sources.

Figure 6 :

 6 Figure 6: γ-ray spectroscopy: estimated 22 Na radionuclide with GMCA, ODL and dGMCA equipped with the robust Fréchet mean. Errors with respect to the input spectrum are displayed in transparent solid lines.

  (a) SNR = 40dB (b) SNR = 80dB

Figure 7 :

 7 Figure 7: γ-ray spectroscopy: evolution of the mixing matrix criterion as a function of the mini-batch size for SNR 40dB and 80dB.

Figure 8 :

 8 Figure 8: Histogram of log 10 of the mixing matrix criterion across mini-batches with the Fréchet mean (left) and robust Fréchet mean (right) and B = 500 and B = 10.

Figure 9 :

 9 Figure 9: Evolution of the mixing matrix criterion as a function of the mini-batches size for a fixed number of batches of 10, 100 and 1000.

This could be easily extended to non-white Gaussian noise by taking into account the noise covariance matrix.

It can possibly vary during the optimization process. It will be kept fixed in this article.

In practice, Nesterov's smoothing parameter was fixed to 0.1, which provides a good balance between speed of convergence and approximation accuracy for this application.

For instance, using the Median Absolute Deviation (MAD) estimator.

Note that in GD, mini-batch optimization can be implemented in a rather natural way as it requires no aggregation step; the computation of the gradient naturally accumulates information from the different mini-batches.

Acknowledgement

This work is supported by the European Community through the grant LENA (ERC StG -contract no. 678282).


Jerome Bobin

Fabio Acero

Adrien Picquenot

METRIC LEARNING FOR SEMI-SUPERVISED SPARSE SOURCE SEPARATION WITH SPECTRAL EXAMPLES

Keywords: blind source separation, sparsity, metric learning, physical constraint

Sparse Blind source separation (sBSS) is an unsupervised matrix factorization problem, which is now a key tool to analyse multispectral data, especially in astrophysics. Being an ill-posed problem, designing priors is crucial. In the present paper, we investigate how the prior knowledge based on examples of physical spectra can be exploited in sBSS, based on the projection onto the barycentric span of these examples. For that purpose, we investigate different metrics to build such projections, and further introduce a novel machine learning approach to build physically relevant reconstruction. Secondly, we show how this can be deployed to design a semi-blind sparse BSS method coined. Preliminary numerical results on realistic astrophysical X-ray images show very promising separation results.

INTRODUCTION

Multi/hyperspectral data are formed of m observations X i , each of which is made of t samples. The standard mixture model defines each observation as a linear combination of n sources S j to which Gaussian noise is added:X = AS + N, where X ∈ R m×t is the data matrix, S ∈ R n×t the source matrix, A ∈ R m×n the mixing matrix and N ∈ R m×t for the noise contribution. BSS aims at recovering both the mixing matrix A and the sources S from the data X only, which is an unsupervised matrix factorization. Since it is ill-posed, it requires additional assumptions about the sources and/or the mixing matrix, such as statistical independence [START_REF] Comon | Handbook of blind source separation[END_REF], nonnegativity [START_REF] Gillis | Accelerated Multiplicative Updates and Hierarchical ALS Algorithms for Nonnegative Matrix Factorization[END_REF] or sparsity [START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF] to only name three. In a very large number of physical applications, examples of spectra (i.e. columns of A) are available that can originate from simulations of parametric models or measurements. The goal of this article is to investigate how these examples can be exploited to constrain the estimated mixing matrix. To the best of our knowledge, this approach has never been investigated so far in the scope of BSS. The proposed approach first This work was supported by the European Research Council -Contract 678282 builds upon constraining the estimated mixing matrix to lie in the barycentric span of the available examples, according to given metric. Several metrics will be studied. Different physical examples will be considered, for which fixed metrics do not yield physically relevant estimates. For that purpose we introduce in Section 2 a machine learning approach that allows learning such relevant metric based on some training set. In Section 3, we show how this example-based constraint can be exploited in the semi-blind BSS setting.

BUILDING MODEL-BASED CONSTRAINTS

Barycentric span and orthogonal projection: In the next, it is assumed that some column a of the mixing matrix can be described by some model M, for which d examples {ϕ e M } e=1,••• ,d are known. Each example is normalized to have unit 1 norm. If the model is parametric, it could be possible to directly estimate the underlying parameters. However, this has two major drawbacks: i) the model can be rather complex to cope with as it can include highly non-linear physical effects as well as instrumental responses, which general leads to highly challenging non-linear optimization problems, and ii) the model may not be known and only some examples of spectra are available. To alleviate these caveats, we propose to constrain a to belong to the barycentric span of the examples according to some metric φ. Formally, this assumption means that there exists some positive weights {λ e } e=1,••• ,d and e λ e = 1 so that a is the solution to:

min a d e=1 λ e φ(a, ϕ e M ), (1)
Several metrics (distances or divergences) can be considered. We will investigate the following three standard metrics: i) the Euclidean distance, which leads to barycenters that are the weighted sum of the examples, ii) the Kullback-Leibler (KL) divergence and iii) the Wasserstein or Earth Mover's distance [START_REF] Villani | Optimal Transport: Old and New[END_REF], which has recently attracted a lot of interest (see [START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised nonlinear dictionary learning[END_REF] and references therein). With the Euclidean and KL metrics, computing barycenters, as defined in Eq. 1, can be tackled using a proximal gradient descent algorithm [START_REF] Parikh | Proximal algorithms[END_REF]. Following [START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised nonlinear dictionary learning[END_REF], Wasserstein barycenters can be efficiently evaluated thanks to an entropic regularization of the original Wasserstein distance and the use of the Sinkhorn iterative algorithm [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF].

In the sequel, enforcing a to belong to some barycentric span requires computing its orthogonal projection onto the barycentric span B φ ({ϕ e M }), which requires finding the weights {λ e } e=1,••• ,d so that:

min {λe} J ({λ e }) = a -b M ({λ e }) 2 2 , (2)
where b M ({λ e }) is the barycenter that is the solution to Eq. 1. The function J ({λ e }) being a smooth, the orthogonal projection P B φ ({ϕ e M }) (a) can be evaluated using a standard gradient descent algorithm with fixed stepsize τ . In practice, evaluating the gradient with respect to the barycentric weights {λ e } e can be tackled using Theano's auto-differentiation function 1 as done for the Wasserstein distance in [START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised nonlinear dictionary learning[END_REF]. A similar strategy also applies to the Euclidean distance and the Kullback-Leibler divergence. In this article, preliminary results have been obtained with fixed gradient stepsize τ (in the range [1e -4, 1e -3]), whose value has been adapted to the metric.

Learning a data-driven metric

The ability of fixed metrics to build physically relevant estimations of spectra will highly depend on the shape of the spectra as well as the evolution of their shape when their underlying physical characteristics evolve. In this article, two distinct meaningful cases will be investigated in the field of X-ray astrophysical imaging: i) emission lines are the instrumental response to a single Dirac at a fixed wavelength, which are generally well described by diracs convolved with a Gaussian-shaped kernel, with a wavelength-dependent width ii) thermal emissions, which lead to observed spectra with several Gaussian-shaped emission lines on top of a continuum. In astrophysical X-ray imaging, the evolution of the continuum yields spectra with different slopes and maximum positions. The amplitude of the emission lines may vary. An example is featured by the purple solid line of Figure 1. In these cases, improving the accuracy of the proposed barycentric constraint requires designing data-driven metrics, which are more adapted to capture the physically meaningful variations of the sought-after spectra. For that purpose, we introduce a new approach that consists in learning metrics that are adapted to build physically relevant barycenters from examples, which makes perfect sense in a large number of applications where training/validation sets can be built based on physical/instrumental models. Let Ω M be some training set made of T samples for a given model M. The proposed metric learning method builds upon the structure of a standard autoencoder [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representationsina deep network with a local denoising criterion[END_REF], which is defined by an encoder Φ that maps the input to the feature space, and a decoder Ψ that maps back to the input space. In contrast to standard autoencoders, metric learning is performed by minimizing the error between the input and the reconstructed projection (according to the Euclidean distance) onto the subspace spanned by the examples i.e. in the learnt/encoded feature space. More precisely, for some encoder Φ, the representation of some input ω ∈ Ω M (resp. some example ϕ e ) is denoted by Φ (ω) (resp. Φ (ϕ e )). One can then define the orthogonal projection of Φ (ω) onto the span of {Φ (ϕ e )} e . To that end, we resort to a simple scheme that further ease the learning process:

• The reconstructed spectra are not enforced to belong to the simplex (i.e. with positive entries and unit 1 norm).

Since the reconstruction error in the input domain is minimized, it should be close to be the case.

• The interpolation in the feature space is not enforced to be a convex combination of the transformed examples. It is substituted by the orthogonal projection onto the subspace spanned by {Φ(ϕ e )} e , which admits a closed-form expression. The weights are then obtained by minimizing the Euclidean distance between the input spectrum and the examples in the feature space:

{ λe (ω)} e = argmin {λe}e Φ(ω) - e λ e Φ(ϕ e ) 2 , (3)
which is a standard least-square problem for fixed Φ. The reconstruction is then defined as b M ( λe (ω)) = Ψ e λe (ω)Φ(ϕ e )

Learning the metric means estimating Φ and Ψ as follows:

min Φ,Ψ ω∈Ω M ω -Ψ e λe (ω)Φ(ϕ e ) 2 2 , (4)
where the weights { λe (ω)} e are defined in Eq. 3. In this article both the encoder and decoder are implemented with multilayer neural networks, such as in standard autoencoders [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representationsina deep network with a local denoising criterion[END_REF].

The number of layers is L, the number of hidden units per layer is fixed to m and activation functions are hyperbolic tangents. The learning stage described in Eq. 4 then boils down to estimating the weight matrices and biases of each layer of both Φ and Ψ. For that purpose, we use the autograd python module2 along with the Adam optimizer.

Illustration In this section, fixed metrics as well as the proposed metric learning approach are evaluated on two distinct cases of spectra that are generally used to describe astrophysical X-ray spectra in the energy range 2 -6keV. These two types of spectra are the Gaussian-shaped emission lines and thermal emissions (see Fig. 1). In these experiments, the number of channels is fixed to m = 50 and d = 2 examples have been selected. For both models 24 samples are available. Amongst them 2 significantly different spectra have been selected at random as examples, 15 have been chosen at random to compose the training set and the remaining 7 form the test set. The same number of training and test samples have been chosen for the two types of spectra. For the fixed metrics (i.e. Euclidean, KL and Wasserstein), the optimization parameters are fixed as follows: i) barycenter evaluation: the step-size is fixed to 10 -3 and the number of iterations to 500, and ii) computation of the orthogonal projection: the gradient stepsize is set to 10 -4 and the number of iterations to 2000. These quantities have been set empirically so that they provide good results for both types of spectra. As for the metric learning method, the number of iterations is fixed to 2.10 4 and the stepsize to 1e -4. These values numerically led to solutions with stable values of the cost function for all tested number of layers and for the various kinds of spectra. The number of layers is set to L = 4 and the number of hidden units per layer to m = 50. In the case of the thermal emission (see Figure 1), none of the fixed metrics are able to capture with accuracy the peak height/position, the slope of the spectrum at high energy and the emission lines. Interestingly, as testified by the reconstruction error (right panel), the metric learning approach leads to a very precise reconstruction of these different structures from only d = 2 examples. More quantitatively, the mean-squared error of the reconstruction from the test set is reported in Table 1. It includes the results obtained for the metric learning method for different numbers of layers. This confirms that learning the metric significantly improves the reconstruction results. Additionally, the accuracy increases with the number of layers up to L = 4; it then collapses for a deeper network. This could very likely be fixed by resorting to standard techniques such as batch normalization or residual network structure, which are known to ease the learning of deep networks. These investigations are left to future work.

SEMI-BLIND SPARSE BSS

The sGMCA algorithm ' In this section, we investigate how the proposed example-based constraint can be combined with sparse BSS to provide more physical and more accurate separation results. To that end, let assume that amongst the n sought-after sources some are perfectly unknown and should be treated blindly and others, say n admit spectra that can be described with n different physical models

{M η } η=1,••• ,n .
The indices of these sources and their corresponding mixing matrix column is denoted by I. In this article, we will only consider cases where a single model applies to a single source; relaxing this limitation is left for future work. Each spectrum A j for j ∈ I is then enforced to belong to the barycentric span of the examples related to model M j as follows:

argmin A,S Λ SW 1 + j =I ι Sm A j + • • • • • • j∈I ι B φ j ({ϕ e M j }) A j + 1 2 X -AS 2 F , (5)
where ι Sm is the characteristic function of the m-dimensional hypersphere (i.e. the corresponding columns of A are enforced to have unit 2 norm), ι B φ j ({ϕ e M j }) is the characteristic function of the barycentric span of the examples {ϕ e Mj } according to the metric φ j and W is the operator that is associated with some sparsifying transform. The problem in Eq. 5 is globally nonconvex but is convex with respect to S. Convexity with respect to A is not ensured since the constraint onto the hypersphere is not rigourously convex and, equipped with the learnt metric, the convexity of the example-based constraint is not guaranteed as well. However, orthogonal projections onto these constraints can be defined and used in an optimization procedure. While several methods can be considered optimization, we opted for a projected alternate least-squares method such as GMCA [START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF], which allows designing algorithms that have a low computational burden, with easy-to-tune parameters. Extending the GMCA algorithm to the proposed semi-blind setting leads to a two-step algorithm semi-blind GMCA (sGMCA), which is detailed in Alg. 1. More precisely, step 1 is composed of two stages: i) the least-squares estimation of the sources, which boils down to applying the left pseudo-inverse of the current estimate of the mixing matrix to the data as in step 1 -a, ii) apply the soft-thresholding operator T with parameter Λ i to each source bf S i . The determination of the values for Λ are described in the next paragraph. The second step copes with the update of the mixing matrix, which is formed of three stages: i) computing the least-squares estimate of A as in step 2 -a, ii) projection each column of A with indices in I onto the barycentric span B φj ({ϕ e Mj }) (described in 2) and iii) normalize the remaining columns to have unit 2 norm.

Parameter setting The initialization is performed using the GMCA algorithm [START_REF] Bobin | Blind Source Separation : the Sparsity Revolution[END_REF]. The number of sources n is assumed to be known and we assume in this article that only a single column of the mixing matrix is described by a single model. Hence, constraining each model to apply to only a single source helps avoiding extra degeneracies. This assumption already allows to tackle a wide range of application scenarios and can be generalized. Then, initialization includes identifying the indices I that associate each model to a single source from the first guess. For each model, this is carried out by finding the column of mixing matrix that provides the minimum Euclidean distance with respect to its projection onto the barycentric span. As emphasized in [START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF], the role of sparsity in BSS is to improve the separation process by increasing the contrast between the sought-after sources, and to increase its robustness with respect to noise. At a given iteration, the error that contaminates a given source estimate is mainly comprised of noise and interferences between the sources, which comes from an imperfect estimation of the mixing matrix. The goal of the thresholding step is then to remove a large part of these two terms. Assuming that they are morphology dense and therefore non-sparse, the standard deviation of this combination of errors can be empirically evaluated using the Median Absolute Deviation estimator (MAD). The regularization parameters Λ can include both a noise-related value based on the MAD and a signal-related term that limits the thresholding bias in the spirit of re-weighted 1 regularization. In this article, the threshold for each source S j will be constant and equal to 3 σ j , where σ j is the MAD-based estimate of the contamination standard deviation (see [START_REF] Kervazo | Heuristics for efficient sparse blind source separation[END_REF] for more details). The maximum number of iterations is fixed to 1000. The stopping criterion is based on the spectral angular distance between consecutive iterates of the mixing matrix. More precisely, for normalized columns of A with unit 2 norm, the algorithm stops when max j arcsin | < A j (k) , A j (k+1) > ≤

, where = 1e -6

NUMERICAL EXPERIMENTS

Experimental setting In this section, comparisons are carried out on realistic simulations of X-ray images provided by the Chandra space telescope 3 . These data correspond to observations of the Cassiopea A supernova remnants [START_REF] Picquenot | Novel method for component separation of extended sources in x-ray astronomy[END_REF].

They are comprised of a linear combination of 3 astrophysical components (displayed in Fig. 3): synchrotron emission, a thermal emission and a redshifted iron (Fe) emission lines (described by the Gaussian-shape emission line) as shown in Fig. 2. These components are representative of typical supernovae remnants in the energy band 5000 -6000 eV (electron-volt) -see [START_REF] Picquenot | Novel method for component separation of extended sources in x-ray astronomy[END_REF]. In the next experiments, the data are composed of 50 observations of size 128 × 128 pixels.

In addition to the sGMCA algorithm, comparisons will be performed with the HALS algorithm [START_REF] Gillis | Accelerated Multiplicative Updates and Hierarchical ALS Algorithms for Nonnegative Matrix Factorization[END_REF], which is a NMF algorithm that minimizes a standard quadratic data fidelity term under positivity constraints. In both the GMCA and sGMCA algorithms, the signal representation W is chosen as the isotropic undecimated wavelet transform [START_REF] Starck | The undecimated wavelet decomposition and its reconstruction[END_REF]. The sGMCA algorithm is tested with two distinct metrics: i) the KL divergence, and ii) the metric learnt with the machine learning approach introduced in Section 2. The metric-based constraint applies to the thermal and line emissions.

In the following comparisons are carried out based on the mixing matrix criterion C A , which is defined as: C A = QA + A -I 1 , where QA + is the pseudo-inverse of the estimated mixing matrix that has been corrected for the scale/permutation indeterminacies with the matrix Q . A is the true mixing matrix. Low values of C A then relates to better separation performances. Each result has been obtained given as the mean of 10 Monte-Carlo simulations with 3 http://chandra.harvard.edu/ different noise realisations.

Results

The contribution of the synchrotron emission is estimated blindly. Furthermore, it is often a dominant emission that hampers the extraction of dimmer emissions. In this paragraph, the robustness of the methods is assessed when the relative flux (ratio with respect to the total 2 norm of the contribution of the other sources) of the synchrotron emission increases from 1 to 20; the SNR is fixed to 15dB. As displayed in Fig. 4, the sGMCA algorithm leads to an improvement of about one order of magnitude, even when the synchrotron emission is highly dominant.

CONCLUSION

We introduce a novel semi-blind source separation method, based on examples. In the proposed approach, semi-blindness is exploited by enforcing the sought-after mixing matrix to have columns that belong to some example-based barycentric span. While fixed metric do not generally yield physically relevant spectra, we further introduce a machine learning approach -inspired by the autoencoders -to learn a metric that is adapted to build accurate and relevant estimations from examples. This constraint is employed in the proposed sGMCA algorithm to tackle sparse semi-blind source separation. Preliminary results based on realistic X-ray astrophysical data, show very promising results, which will be further detailed during the conference. Initialization based on GMCA 3:

L =
while Convergence is not reached do 4:

(1-a) Least-square estimation: S = A (k) + X 5:

(1-b) Thresholding: S (k+1) = T Λ (S W) W T 6:

(2-a) Least-square estimation: A = XS (k+1) + 7:

(2-b) Projection onto the barycentric span:

8 :b

 8 ∀j ∈ I, A j (k+1) = P B φ j ({ϕ e M ({λ e }) = argmin a d e=1 λ e φ(a, ϕ e M ) 12: end procedure

Fig. 1 :

 1 Fig. 1: Left: reconstructed thermal emission spectra of the KL, Wasserstein, Euclidean and learnt metrics. Right: estimation errors.

Fig. 2 :

 2 Fig. 2: Left: input and estimated spectra with sGMCA. Right: Estimation errors for HALS and dGMCA with the KL and the learnt metric.

Fig. 3 :

 3 Fig. 3: Left: input synchrotron source. Middle: input iron (with emission line spectrum). Right: observation at 5548eV.

Fig. 4 :

 4 Fig. 4: Evolution of the mixing matrix criterion as a function of the synchrotron relative amplitude.

Table 1 :

 1 MSE in dB of the reconstructed spectra from different metrics.

		3	4	5	Eucl. Wass. KL
	G-shaped	6.7	39.3 4.6	-0.5	13.5	25.4
	Thermal	40.5	47.2 21.8 13.5	10.6	13.7
	Algorithm 1 semi-blind GMCA
	1: procedure
	2:

https://github.com/Theano/Theano

https://pypi.org/project/autograd/


Abstract--The detection of single line-to-ground (SLG) fault in compensated distribution networks can be hampered by important capacitive currents. The way of dealing this problem depends on particular regime where data is acquired. In transient regime the relays diagnosis based on analysis of charging components can deceive in massive presence of discharging currents. The latter however will disclose data for directional function, if sampling frequency is high enough for acquisition in wave propagation area.

The fault detection can also be out of reach of steady state methods which are based on detection of active component in faulty currents. Efficiency of these methods vanishes when strong capacitive currents squeeze the phase shift between faulty and sound zero sequence currents below actual acquisition limits. In such situation however it is still possible to get the faulty feeder indication from the way the residual currents cross the zero level soon after inception.

Index Terms--Fault diagnosis, power distribution protection, compensated distribution system

I. INTRODUCTION

N case of low ohmic SLG faults in compensated distribution networks the directional function can be assured by transient relays which compare polarities of zero sequence charging components at main frequency, either of voltage versus currents on each feeder [START_REF] Nikander | Applications of transient phenomena during earth faults in electricity distribution networks[END_REF] or of their products on different feeders [START_REF] Coemans | Using the EMTP and the Omicron to design a transients based digital ground-fault relay for isolated or compensated networks[END_REF]. Some difficulties can arise however in presence of violent discharging transients. They grow with inception angle θ, which determine pre-fault values of the faulty line voltage. The discharging currents are absent at θ=0° but then they reach the charging currents' level when the inception angle is several degrees (Fig. 1), and grow rapidly. The most favorable conditions for generation of important transients are with low resistance fault when it occurs at inception angle 90° in a capacitive network, be it a cabled system or a mixed one, composed of cables and lines.

In such conditions we can expect the extraction of main frequency component out of the charging currents to be more delicate an operation in cables than in lines. The charging This work was supported in part by Tsing Hua University in Taiwan and by Schneider Electric, France Wan-Ying Huang is with Supelec, Gif sur Yvette, France (phone:33 1 69 85 15 32; fax: 33 1 69 85 15 39; e-mail: wanying.huang@supelec.fr) Robert Kaczmarek is with Supelec, Gif sur Yvette, France (e-mail: robert.kaczmarek@supelec.fr) components in overhead lines [START_REF] Lehtonen | Fault management in electrical distribution systems[END_REF] were reported to reach amplitudes 10-15 times the rated frequency amplitude whereas the discharging currents were estimated as several per cent of charging components. However, when replacing lines with cables, all the other conditions unchanged, the discharging currents will be more important than the charging ones, with frequency span between them diminishing.

A relevant example comparing amplitudes and frequencies of the main discharging and charging currents (Fig. 2) presents the amplitudes ratio A dis /A ch rising from 0.25 (lines) to 1.5 (cables) and the frequency ratio f dis /f ch diminishing from 11 to 6. The reason of these tendencies is with specific values of zero sequence capacitances in cables and lines.

Obviously, it is simpler a task to isolate a paramount and somewhat isolated component. The consequence for the transient relays, when applied in capacitive systems, can be an uncertain choice of window and trigger for acquisition and heavy standards on extraction of the relevant charging components.

SLG Fault Detection in Presence of Strong Capacitive Currents in Compensated Networks

Wan-Ying Huang and Robert Kaczmarek I Fig. 2. The faulty feeder zero sequence currents in cables and lines of the same length and power transfer, with 1Ω fault resistance and 90° inception angle. The discharging components in cables are of higher amplitudes than the charging components and the frequency ratio are different in cables and in lines.

The disturbing presence of the discharge components can however turn into an opportunity to make a correct directional decision. This opportunity is offered by rigorous waveform disposition in initial propagation zone, usually unexplored in distribution systems. However, as the high frequency acquisition procedures become simpler and cheaper, we are tempted with evaluation of possibilities of using the travelling wave regime to get directional function.

Strong capacitive currents are troublesome also in case of resistive faults, where relays discriminate faulty and sound feeders upon presence of active current component in residual currents. This can be achieved either by comparing signs of projections of zero sequence currents [START_REF] Griffel | Anew deal for safety and quality on MV networks[END_REF][START_REF] Welfonder | Location strategies ans evaluation of detection algorithms for earth faults in compensated[END_REF] or looking for phase advance [6,7] of the faulty zero sequence current over the sound ones. However, in steady state the strong capacitive currents diminish this phase advance, with possible inhibition of the discrimination capacities of relays.

Then the way to re-establish these capacities is to exploit data recorded in transient regime, where the apparent phase difference is more important than in the steady state.

The procedures which follow were modeled and simulated in EMTP using frequency dependant parameters.

II. DIRECTIONAL FUNCTION OF DISCHARGING CURRENTS

We consider a radial network (Fig. 3) with no discontinuities in feeders' impedances. The network is supplied through a transformer with secondary winding grounded through Petersen coil. An SLG fault is modeled as a resistance R f . Fig. 3. A phase-to-ground fault in a radial network, all the sound feeders aggregated into one At fault occurrence an initial voltage wave annulling the phase voltage travels along the faulty phase of the faulty feeder with negative amplitude, accompanied by current wave also of negative amplitude. Their shape will be modeled by multiple refractions and reflections from busbar, fault and loads. The faulty phase current arriving on busbar reflects and the resulting current is equally distributed among faulty phases of all the sound feeders. These faulty phase currents impose their waveform upon residuals in each feeder. Consequently, initial currents on the sound feeders will be measured with the same polarity and opposed to the polarity on the faulty feeder.

On each sound feeder the initial polarity regime will be over with first busbar reflection of wave getting back from loads. Its overall travel time is 2l/v, with l -the feeder's length and v -the wave velocity. This is the time interval where the current on the faulty phase of any sound feeder keeps its initial polarity unchanged. The shortest distance l ss is with the shortest sound feeder and thus we get the duration of the initial polarity v l ss ip / 2 = δ , where v is the maximal modal velocity.

The parameter

ip ip f δ / 1 =
determines the minimal sampling frequency necessary to get one point in the zone of a rigorous polarity disposition of residual currents. For example, with the shortest length of a cable feeder being 2km this frequency can be of 20 -30 kHz. For practical reasons the frequency of about 100kHz is to be reckoned with.

In cables the propagation phenomena take place both in cores and sheaths. The fault we refer to is a core-to-sheath (and, eventually, -to-ground) piercing, with sheaths grounded on supply side or on both sides. In each of these cases the traveling waves are very similar and the initial polarity zone clearly exposed (Fig. 4). The initial polarity time interval does not depend neither on fault position, nor fault resistance value. If scrutinized on all feeders it can point out from busbar to fault in systems with laterals [8]. In order to start tracing of the faulty feeder, the initial polarities of residual currents in all busbar connected feeders have to be compared. The beginning of the faulty chain is pointed with a unique sign, called "witness sign", being different from polarity of all the other feeders. Then we follow current sensors on feeders with the same "witness" polarity. The fault is at the end of the chain.

This will be illustrated in the case of an SLG fault on the feeder d in a five feeders network (Fig. 5).

The rising current profiles have been recorded first on feeders d and e (Fig. 6). When the traveling waves arrive to busbar, the sign of the current measured on the feeder c is different than those measured on the feeders a et b. The feeder c disclosing the "witness sign", it points out to the feeder d as the fault location. All methods based on analysis of traveling waves are highly sensitive to impedance mismatches, what results in severe conditions on their application in distribution networks. Feeders with single tee joints can be analyzed relatively easily, unless the initial polarity zone is too short to be detected. On the contrary, multiple tee joints and joints between cables and lines in mixed systems make the detection problem inextricable in terms of traveling waves' analysis.

III. THE PHASE DIFFERENCE OF THE ZERO SEQUENCE CURRENTS

A. Principle

The steady state SLG fault regime can be analyzed on equivalent residual circuit (Fig. 7), where V'' is the voltage over a SLG fault emplacement in absence of the fault, I 0_fzero sequence current on the faulty feeder, ΣΙ 0_s -the sum of zero sequence currents on all the sound feeders and I N is the neutral point current composed of a resistive and an inductive components.

During permanent fault regime, an active current component is present in zero sequence current I 0_f on the faulty feeder. The resulting phase difference between the faulty and the sound zero sequence currents I 0_s (Fig. 8) is the basis of traditional wattmetric method of detection.

In low capacitive lines the phase advance can be almost 90°, because under the effect of compensation we have (Fig. 7):

( )

1 _ _ 0 _ 0 RN f c s N f I I I I I - - = - - =

∑

where the faulty residual I 0_f is dominated by its active component I RN . It is then easy to take direction decision. In cables however the faulty feeder capacitive current I c_f dominates the composition of I 0_f and diminishes readability of the phase advance, particularly with fault on a long feeder or in case of system over tuning.

We have simulated a three feeders' cable system of capacitive currents 15A+10A+5A, the fault installed on the 15A feeder, with 100% tuning and the fault resistance values ranging from 1Ω to 1000Ω. In all cases the actual phase advance in steady state was less than one sample step at 600Hz sampling frequency, this frequency being used in certain relays.

Fortunately, when the phase advance in steady state becomes undetectable, we can look upon an analogous parameter in transient regime, where it can be much larger. This is a consequence of the way the transient regime develops, beginning just after the fault inception with phase opposition between faulty and sound current residuals and finishing in steady state with a slight phase advance of the faulty residual over all the sound ones. This development is correlated with evolution of the neutral point current I N smoothly growing from zero to its permanent value. During the first millisecond after fault inception it can grow very slowly, particularly with resistive faults, because of high values of the neutral point elements L N and R N .

On the contrary, the feeders' line-to-ground capacitors C f and C s charge and discharge vigorously. During a short time interval the neutral point current I N is negligible comparing to capacitive zero sequence currents.

The latter being under the same charging conditions, the faulty zero sequence current is initially in phase opposition to the sound feeders zero sequence currents, see (2):

( ) ( ) 2 _ 0 _ 0 ∑ - = s f I I
and proceeds toward zero level (Fig. 9) with different polarities.

The identification of faulty feeder operates then with aid of following algorithm.

B. Algorithm

We detect the slopes of filtered residuals at their first zero crossing after the fault inception. If all but one witness the same slope sign, then we can declare the latter as the faulty one without even controlling its zero crossing:

=         =         =         - = = - L L
This is a one shoot procedure, without possibility of verification. On the other hand it is a conclusive test, as the matter goes about unambiguous identification of slopes' signs.

IV. CONCLUSION

We think that strong capacitive currents, sometimes generating unfavorable conditions for protection relays in compensated systems, can be exploited as carriers of relevant information.

Whenever extremely rapid information is required, we can find it when treating the discharging currents as useful, at a price of higher sampling frequency. After fault inception we dispose of several tens of microseconds to get the data in wave propagation area. This can make sense in simple distribution systems, with single ramification per feeder and homogenous line impedance.

We can also make useful the presence of strong capacitive currents in permanent fault regime, where these currents occult detection of faulty feeder. In such cases the diagnosis based on phase advance of the faulty residual current over the sound residuals is better assured when tracing the corresponding apparent phase advance in transient regime.

Both propositions need the current data be centralized, they work on one or only few data points in a very short time span. These drawbacks are price for matching consequences of cable proliferation in resonant grounded distribution networks.

Fig. 1

 1 Fig. 1 Initial faulty feeder zero sequence currents in cables. The discharging currents of higher frequency are superimposed on the charging ones of lower frequency. Upper figure: θ=0°, no discharging components. Bottom: θ=5°, the discharging current's amplitude reach the level of charging components.

Fig. 4 .

 4 Fig. 4. The initial polarity zone in three feeders' network, zero sequence currents. The index s is for "sound" and f-for "faulty" feeder

Fig. 5 .

 5 Fig. 5. A phase-to-ground fault F on the feeder "d" in a network with laterals

Fig. 7 . 3 Fig. 8 .

 738 Fig. 7. An equivalent residual circuit of the faulty network on Fig. 3

Fig. 9

 9 Fig. 9 Faulty and sound current residuals in case of resistive fault, non filtered curves (Rf =1kΩ, θ=90°)

V. ACKNOWLEDGMENT

The authors are grateful to Mr Volker Leitloff for his expert comments on initial draft of this paper.

Robert KACZMAREK is professor in Ecole Supérieure d'Electricité (Supelec) in France. His research interest is with power networks (fault detection, network protection, parameters identification), magnetism (domain observation, loss prediction and measurements in soft and hard magnetic materials) and design of electric machines.