Je souhaiterais tout d'abord remercier les membres du jury, ainsi que mes deux rapporteurs Nicolas Gillis et Nicolas Dobigeon, pour leur questions pertinentes et leur lecture attentive de ce manuscrit.

Mes remerciements vont ensuite à mon directeur de thèse, Jérôme Bobin, pour son aide précieuse et efficace, ainsi que pour les nombreuses discussions qui en ont découlées, me permettant notamment de mieux appréhender le monde de la recherche. Pour ses conseils avisés, je suis également reconnaissant à Florent Sureau, qui m'a permis d'approfondir de manière SURE certains aspects de ce doctorat.

J'ai également beaucoup apprécié le temps passé au sein du groupe CosmoStat. Je remercie son fondateur, Jean-Luc Starck, ainsi que tous ses membres. Je pense plus particulièrement à ceux avec qui j'ai partagé le bureau, avec (par ordre d'arrivée) : Cécile, Joana, Imane et Jiaxin. Parmi les membres de l'équipe, je souhaiterais également particulièrement remercier les deux stagiaires que j'ai eu le plaisir de coencadrer, à savoir Tobias Liaudat et Yarui Zhang, pour leur contribution ainsi que l'apport pédagogique. Enfin, merci à tous ceux avec qui j'ai partagé encore davantage de bons moments : les membres du Cheeseday et les autres doctorants parmi autres.

J'aimerais également remercier les nombreux amis externes au laboratoire qui m'ont soutenu durant ces trois ans : que ce soit ceux d'école primaire, ceux du collège, du lycée, des classes préparatoires ou d'école d'ingénieur, leur amitié qui a su perdurer depuis compte beaucoup à mes yeux. Enfin, j'aimerais dédier ce travail à ma famille et notamment à mes parents. Je pense que le contexte personnel qu'ils m'ont fourni tout du long de mon enfance a pour beaucoup joué dans la possibilité d'arriver jusqu'à une soutenance de doctorat.

« La science a-t-elle promis le bonheur ? Je ne le crois pas. Elle a promis la vérité, et la question est de savoir si l'on fera jamais du bonheur avec la vérité »

Notations and abbreviations

General notations

-u : a scalar. u 1..n is used as a shortcut for writing the whole set of scalars u 1 , u 2 , ..., u n ; u : a vector ;

-U : a matrix. Vectors are considered as a matrix with only one column ;

-U : a set ;

-U i : i th line of matrix U (when U = S corresponds to the sources, also called an observation). The notation is extended to subsets of lines : if J is a finite subset of [1,n], U J denotes the lines indexed by J ; -U j : j th column of matrix U (when U = S corresponds to the sources, also called a sample). The notation is extended to subsets of columns : if J is a subset of [1,t], U J denotes the columns indexed by J ; -U ij : (i,j) th entry of matrix U (also called coefficient) ; -U T : transpose of U ;

-U † : Moore-Penrose pseudo-inverse of U ;

-U * : true underlying matrix (to be estimated) ;

-Ũ : in proximal algorithms, estimate of a variable U * before the application of the proximal operator ;

-Û : estimate of U * by an algorithm ;

-Û(l) : in iterative algorithms, denotes the estimation of U * at iteration l ; Û(1..l) is used as a shorthand for the set of all the estimates until iteration l.

-U[J] : used to denote an estimation of U * using the mini-batch indexed by J. Note that we do not use a hat here, as it is used for the estimate obtained from the aggregation of the various U[J] ; f (.) function with scalar output ; f (.) function with matrix output ;

-1 n×t : matrix of size n × t filled with ones ;

-Id : identity matrix ; -J C complementary set of J ;

-Diag(λ 1 , ..., λ n ) : diagonal matrix having as diagonal elements λ 1 , ..., λ n .

Operators and norms

u p , p ∈ R + : p (quasi)-norm of u, that is for u ∈ R u , u p = ( u i=1 u p i )

Contents

-U p matrix norm of U induced by the p -norm on vectors ;

-U ∞ : maximum absolute value of the coefficients of U. The notation is the same for vectors : u ∞ ; -: Hadamart product (e.g. elementwise product of two matrices) ;

-U 0 means that all the coefficients of U are non-negative : U i,j ≥ 0 ; -< .|. > : scalar product ;

-∆ f (.) : gradient of f ;

-prox f (u) : proximal operator of f in u. See Appendix A for definition and special case of soft-thresholding S λ (u) and projection on set U , Π U (u) ;

-M f (.) : Moreau envelope of f ;

ι U (.) is the characteristic function of a set U :

∀u ∈ R, ι U (u) = 0 if u ∈ U ∞ otherwise
The notation is extended for matrices m : number of observations ;

U ∈ R m×n : ι U (U) = m i=1 n j=1 ι U (U ij
t : number of samples ;

-T : number of samples in a transformed domain ;

σ : a standard deviation ; Contents 3 h : differentiable term in a cost function (most of the time, h is also multiconvex) ;

h (i) : convex function corresponding to the multi-convex function h with all but one block fixed ;

-J (.), G(.) and J (i) : functions used to enforce constraints on A and S respectively, while J (i) is used for the more general case of more than two A and S matrices. Examples of constraints : positive orthant

K + = {S ∈ R n×t ; ∀j ∈ [1, n], k ∈ [1, t], S k j ≥ 0} and oblique constraint O = {A ∈ R m×n ; ∀j ∈ [1, n], A j 2 2 = 1} ; -K : number of constraints J (i) ;
-Φ or Φ S : a sparsifying transform (size T × t) ; -R S (size n × T ) control the trade-off between the data fidelity and the sparsity terms. It can be decomposed into R S = Λ S G where Λ S (n × n) is a diagonal matrix of the regularization parameters λ 1 , λ 2 , ..., λ n and G (n × T ) is a matrix used to introduce individual penalization coefficients in the context of reweighted 1 ;

-M S : GMCA regularization parameters. M S (l) = Diag(µ

(l) 1 , µ (l)
2 , ..., µ

n )1 n×t ;

γ, δ : parameters in (0, 1) used for the step size ;

η : gradient step size, or learning rate in the machine learning terminology ;

-L : Lispchitz constant of a gradient ;

ω i , ∀i ∈ [1, B] : weights : ω i ≥ 0 and B i=1 ω i = 1. The ω i can be written within a vector w : w i = ω i ; κ : constant used as a multiplicative factor of the mad. Usually, κ = 3 ; l : in iterative algorithms, number of the current iteration (l f is specifically used for the Fréchetalgorithm) ;

-L : number of iterations of an iterative algorithm (L f is specifically used for the Fréchetalgorithm) ; Chapitre I

Résumé

La séparation aveugle de sources (BSS 1 - [Comon & Jutten 2010]) est une méthode de premier plan pour apprendre des décompositions physiques de données multi-valuées. Celle-ci a fait ses preuves dans de nombreux domaines, tels que par exemple le traitement du signal audio [START_REF] Vincent | [END_REF], Vincent et al. 2011, Ozerov & Févotte 2010, Duong et al. 2010, Févotte et al. 2009], le biomédical [START_REF]Removing electroencephalographic artifacts by blind source separation[END_REF], Negro et al. 2016, Poh et al. 2010] et l'astrophysique [Bobin et al. 2014].

Les jeux de données utilisés en BSS sont obtenus à partir de mélanges de signaux élémentaires appelés sources (cf. Fig. I.1 pour un exemple illustratif tiré de données astrophysiques). De manière générale, les observations, regroupées en tant que lignes d'une matrice X (de taille m × t), peuvent donc s'écrire comme : X = f * (S * ) + N où les sources sont les lignes de la matrice S * (de taille n×t), qui sont mélangées par la fonction f . La matrice N correspond aux imperfections du modèle et est appelée bruit. Une instance spécifique de ce type de problèmes est le cas où le mélange est linéaire : X = A * S * + N (I.1) Par exemple, les sources peuvent être supposées statistiquement indépendantes, ou avoir des coefficients positifs. Dans ce manuscrit, nous nous intéressons spécifiquement au cas où les sources sont parcimonieuses, puisque cette approche a permis d'obtenir d'excellents résultats lors de la dernière décennie [START_REF] Zibulevsky | [END_REF], Bobin et al. 2007, Bobin et al. 2015]. Pour résumer très succinctement, l'hypothèse de parcimonie présuppose que les sources comportent un grand nombre de coefficients nuls (potentiellement dans un espace transformé, par exemple le domaine de Fourier pour donner un cas simple).

(a) X 1 (b) X 5 (c) X 9 (d) X 30
En tant que tel, le problème de BSS parcimonieuse peut s'écrire comme un problème d'optimisation multi-convexe de la forme (nous prendrons ici le cas le plus simple) :

argmin A∈R m×n ,S∈R n×t 1 2 X -AS 2 F + R S S 1 + ι {∀i∈[1,n]; A i 2 2 =1} (A) (I.2)
Dans lequel :

-Le terme 1 2 X -AS2 F , avec . F la norme de Forbenius, promeut une reconstruction fidèle des données.

-Le terme R S S 1 encourage la parcimonie des sources, où est utilisé pour le produit d'Hadamard. La matrice de paramètres R S (de taille n × t) permet de contrôler le compromis effectué entre les termes d'attache aux données et de parcimonie.

-Pour éviter des solutions A and S dégénérées dans lesquelles A F → ∞ et S F → 0 à cause du terme de parcimonie, un dernier terme appelé contrainte oblique est introduit. Il impose que toutes les colonnes de A se trouvent sur l'hypersphère unité 2 . La fonction charactéristique est notée ι.

Les algorithmes classiques de BSS parcimonieuse s'attachent donc généralement à minimiser le problème (I.2), soit de manière exacte 2 (BCD [Tseng 2001] ou PALM [START_REF] Bolte | [END_REF]), soit de manière approximative mais en introduisant des heuristiques3 permettant d'augmenter empiriquement la robustesse grâce à un choix automatique des paramètres de régularisation R S (GMCA [Bobin et al. 2007]).

Toutefois, en dépit de leurs nombreux succès, la plupart des méthodes ont été utilisées sur des problèmes de petites tailles. Par conséquent, le déluge de données actuel représente un important défi pour les méthodes de BSS actuelles. En astrophysique par exemple, de nouveaux instruments tels que le Square Kilometer Array (SKA) seront à même de fournir des volumes de données colossaux. Dans le contexte de la BSS, et pour donner un ordre d'idée des tailles à considérer, des données comprenant jusqu'à t = 10 9 échantillons et m = 104 observations pourraient faire leur apparition. De manière similaire, en spectroscopie les mélanges peuvent comporter plusieurs dizaines de sources. L'objectif de ce doctorat est donc de proposer de nouvelles méthodes de séparation aveugle de sources parcimonieuses permettant de traiter des problèmes grande échelle.

Plus précisément, les travaux s'articulent autour de quatre problématiques majeures :

1 -L'introduction d'une méthode de choix automatique des paramètres de régularisation R S du problème (I.2) lors de sa minimisation par l'algorithme PALM [START_REF] Bolte | [END_REF]. Cette problématique est primordiale, puisque qu'un tel choix est généralement effectué en testant différentes valeurs de R S et en regardant la qualité de la solution obtenue par PALM pour chaque jeu de paramètres (approche de type grid-search en Anglais). Néanmoins, dans le cas grande échelle, les temps de calculs prohibitifs induits par cette approche sont souvent rédhibitoires, nécessitant l'introduction d'une méthode automatique permettant de ne lancer qu'une unique fois l'algorithme ;

2 -L'introduction d'une méthode permettant de traiter un nombre de sources n important. Le défi présente ici deux facettes, puisque les algorithmes classiques de BSS, lorsqu'ils sont confrontés à un nombre de sources croissant i) nécessitent des temps de calcul accrus ; ii) voient la qualité des solutions qu'ils proposent décliner ;

Après une introduction plus détaillée du problème et une présentation des stratégies d'optimisation utilisées dans les chapitres II et III respectivement, le chapitre IV s'attaque à la résolution du problème 1). Plus spécifiquement, la première partie consiste à déterminer les difficultés induites par l'utilisation de PALM dans le cadre de la BSS parcimonieuse. Une étude empirique nous permet d'avancer que dans ce contexte, utiliser PALM avec une approche de type grid-search pour la recherche de bons 4 paramètres de régularisation R S souffre d'une faible efficacité et versatilité. Pour résumer succinctement, une faible efficacité signifie qu'il est délicat, pour une expérience donnée, de trouver de bons paramètres de régularisation R S . Une faible versatilité implique qu'il est difficile de généraliser un bon choix de R S d'un jeu de données à un autre. Ainsi, faibles efficacité et versatilité rendent le choix des paramètres de régularisation très complexes lors de l'utilisation de PALM pour des problèmes grande échelle 5 . Cependant, l'étude proposée montre qu'une fois de bons paramètres R S trouvés, l'estimation des matrices A * et S * par PALM peut être de très bonne qualité, surpassant celle d'algorithmes heuristiques tels que GMCA. C'est pourquoi nous proposons, dans une deuxième partie, une méthode permettant de contourner les difficultés précédentes. Plus spécifiquement, une approche en deux étapes est introduite, comprenant :

-Une étape d'initialisation se basant sur l'algorithme GMCA : celui-ci, bien qu'approximatif car basé sur des moindres carrés alternés projetés (pALS), propose une méthode automatique de choix de R S , qui donne en pratique une première estimation décente  et Ŝ ;

-Une étape de raffinement se basant sur PALM. L'initialisation correcte fournie par GMCA permet de déterminer un bon jeu de paramètres R S , utilisé dans cette deuxième étape. PALM permet alors potentiellement d'améliorer la solution fournie par GMCA, tout en fournissant des garanties mathématiques, telle la convergence de l'algorithme.

La qualité de l'approche est démontrée grâce à une expérience de BSS sur des données réalistes d'astrophysique. Une discussion est également proposée quant aux limitations de la méthode.

4. Le terme "bon" est à comprendre au sens donnant des estimations  et Ŝ proches des vrais facteurs physiques A * et S * . Cette notion est donc distincte de simplement trouver un minimum local du problème (I.2) : il s'agit ici de trouver un minimum particulier correspondant à des facteurs ayant un sens physique.

5. Et ce d'autant plus que PALM souffre en plus d'une fiabilité limitée, signifiant sa sensibilité à l'initialisation : autrement dit, il faut également potentiellement relancer plusieurs fois l'algorithme avec des initialisations différentes jusqu'à l'obtention de résultats corrects.

Dans le chapitre V, nous nous attaquons au problème 2). Nous montrons sur une expérience simple la difficulté de traiter des cas de BSS comprenant un grand nombre de sources n. Pour traiter ce problème, la solution proposée se concentre sur la stratégie d'optimisation. Plus spécifiquement, nous utilisons une méthode de minimisation par blocs de coordonnées, ce qui est couramment utilisé en BSS pour contourner l'aspect non-convexe de la fonction de coût (I.2) en utilisant sa structure multi-convexe. Les approches classiques peuvent être catégorisées en deux familles :

-Les méthodes de déflations ou hiérarchiques : celles-ci utilisent des blocs de taille 1 : à chaque itération, une seule source est mise à jour, ainsi que la colonne de A correspondante ;

-Les méthodes utilisant l'intégralité des matrices A et S. Dans ces méthodes, telles le GMCA usuel, seulement deux blocs sont utilisés, le premier correspondant à la matrice de mélange et le deuxième à la matrice source.

Par opposition, la méthode proposée, appelée block-GMCA (bGMCA), introduit des blocs de tailles intermédiaires r au sein de l'algorithme GMCA 6 . Ainsi, alors que les approches par déflation ou hiérarchique correspondaient au cas r = 1, et les méthodes utilisant l'intégralité des matrices au cas r = n, bGMCA utilise des tailles r ∈ [1, n].

Outre un gain substantiel en termes de temps de calcul, il est montré expérimentalement que l'algorithme permet une forte amélioration de la qualité de la séparation pour des tailles de blocs modérées. Dans des cas simples, bGMCA permet aussi de retrouver de manière quasi-exacte (aux incertitudes numériques près) les facteurs A * et S * . L'explication avancée pour ces résultats est qu'il existe un compromis dans la taille des blocs : pour r proche de 1, l'algorithme souffre d'erreurs de propagation entre les itérations et le fait d'utiliser des blocs crée des erreurs, assimilables à un bruit supplémentaire. Pour r = n, le problème de séparation devient plus complexe. Prendre des valeurs de r intermédiaires semble donc permettre de limiter les erreurs de propagation, tout en bénéficiant d'un problème plus simple.

Pour conclure ce chapitre et montrer que bGMCA peut être utilisé dans le cadre de la BSS avec des a priori plus complexes que la seule parcimonie, un problème réaliste de LC / 1 H NMR est proposé, dans lequel la positivité des sources et de la matrice de mélange est de surcroit imposée.

Dans le chapitre VI, le problème 3) est pris en considération. Plus spécifiquement, bGMCA permettait de traiter des problèmes avec un grand nombre de sources mais nécessitait l'utilisation de l'ensemble de la matrice de données X à chaque itération. L'algorithme proposé ici, nommé distributed-GMCA (dGMCA), élimine cette 6. Plus exactement, des blocs sont introduits dans l'algorithme en deux étapes décrit précédemment.

limitation. Le principe général consiste en l'introduction de B mini-batchs [START_REF] Xing | [END_REF]] dans le schéma de minimisation par projections alternées des moindres carrés, et plus spécifiquement dans GMCA, ce qui permet de continuer à bénéficier de la grande robustesse de cet algorithme, ainsi que de son choix simple de R S 7 . Cependant, l'utilisation du pALS implique qu'à chaque itération, un nombre B d'estimations de la matrice de mélange A * soit calculé, soit une par mini-batch. Ces estimations ne sont pas toutes d'égales qualités, puisque réalisées sur des parties différentes, de petite taille, des données X, qui peuvent être plus ou moins simples à démélanger. Une question naturelle est donc comment les agréger, pour obtenir à chaque itération une unique bonne estimation de A * . Une méthode naturelle serait de réaliser l'agrégation en effectuant une simple moyenne euclidienne. Cependant, l'estimée de A * doit également respecter la contrainte oblique présente dans (I.2), ce qui n'est pas garanti par l'utilisation de la moyenne euclidienne. Nous proposons donc plutôt d'utiliser une moyenne de Fréchet, c'est à dire une moyenne sur l'hypersphère. Pour robustifier l'approche, une deuxième version de l'agrégation est proposée, similaire à une médiane sur l'hypersphère. De plus, une pondération est utilisée, permettant de prendre en compte la qualité de l'estimation des sources dans chaque mini-batch. Les expériences réalisées permettent de démontrer le gain de temps (en plus du gain en mémoire, puisque X est potentiellement découpé sur différents noeuds d'un cluster) obtenu par la méthode. L'étude de la qualité de la séparation permet par ailleurs de distinguer deux régimes :

-Dans le cas de sources modérément parcimonieuses, la séparation obtenue lors de l'utilisation de mini-batch est quasiment identique à celle en utilisant l'ensemble des données, si tant est qu'une taille de mini-batch raisonnable soit utilisée. Introduire des mini-batchs permet donc, à coût quasiment nul en termes de qualité de séparation, un gain en temps de calcul important ;

-Dans le cas de sources très parcimonieuses, les résultats sont plus surprenants : pour des mini-batchs de faibles tailles, la séparation est même meilleure que lorsque l'ensemble des données est utilisée. L'explication proposée se base sur des liens avec certains travaux récents sur la descente de gradient stochastique (SGD). Notre hypothèse est que, similairement à la SGD, l'introduction de mini-batch favorise certains minima du paysage d'optimisation de la fonction de coût (I.2) qui généralisent bien. Dans le contexte de la BSS, de tels minima correspondraient à des solutions peu sensibles à une réalisation donnée des sources.

Ces résultats sont confirmés sur une expérience réaliste en spectroscopie gamma.

Le chapitre VII est une extension des travaux précédents, et s'attaque au pro-7. En effet, bien que le chapitre IV ait permis de trouver une méthode pour fixer les paramètres de régularisation dans PALM, il est à noter que l'approche nécessite au moins quelques itérations de GMCA pour initialiser l'étape de raffinement. Par conséquent, faire fonctionner le pALS dans le cas grande échelle reste nécessaire. blème 4). Peu d'études se sont intéressées à la séparation aveugle de sources nonlinéaire. Ceci est probablement dû à des indéterminations beaucoup plus importantes que dans le cas linéaire. Notamment, l'indépendance des sources n'est par exemple plus suffisante pour garantir leur séparation. Qui plus est, même dans le cas où les sources seraient séparées, le reconstruction n'est pas garantie puisque celles-ci sont retrouvées à une fonction non-linéaire h près, qui ne les remélange pas. L'algorithme proposé, nommé StackedAMCA, bénéficie d'une interprétation géométrique du problème de BSS non-linéaire : graphiquement, chaque source est transformée en une variété de dimension 1 (1D). Nous proposons d'estimer les non-linéarités en ajustant un modèle linéaire par morceaux à ces variétés grâce à un algorithme itératif. A chaque itération, un nouveau morceau est estimé grâce à un algorithme de BSS parcimonieuse robuste aux non-linéarités de grandes amplitudes [Bobin et al. 2015]. Pour passer d'un morceau à l'autre, et ainsi préparer l'itération suivante de l'algorithme, un résidu est utilisé, correspondant aux données X desquelles sont soustraites les contributions des modèles linéaires déjà estimés. Pour limiter les erreurs de propagation, à chaque itération l'algorithme repart des données originales X et déploie les variétés en appliquant l'ensemble des inverses des modèles linéaires calculés auparavant. L'algorithme admet de plus une interprétation en termes de réseau de neurones, qui est détaillée. La partie expérimentale comporte plusieurs expériences : un mélange linéaire par morceaux, permettant ainsi d'étudier les mécanismes de StackedAMCA dans un cas où le processus de démélange correspond exactement à celui du mélange ; un mélange compliqué dans lesquelles sont présentes de nombreuses sources : Stacke-dAMCA obtient une meilleure qualité de séparation que les autres méthodes de l'état de l'art. De manière intéressante, la méthode est aussi ici capable de reconstruire les sources, malgré l'indétermination par la fonction h : dit autrement, la structure de l'algorithme crée une régularisation implicite qui permet la reconstruction. Nous proposons une dernière expérience dans laquelle la reconstruction n'est plus garantie, dans l'objectif de montrer que l'algorithme sépare encore toutefois bien les sources. Enfin, dans une dernière partie, les hypothèses requises pour StackedAMCA sont étudiées, et en particulier un sous-ensemble de mélanges pour lesquels l'algorithme est supposé être capable de reconstruire les sources est caractérisé.

Chapitre II

Sparse Modelling and Large-Scale BSS

A Multi-valued data analysis and BSS

The overwhelming quantities of data collected everyday have lead to the current so-called deluge of data. Such a phenomenon finds its roots in many origins : one can for instance cite the development of social media and search engines, or new industrial equipment, cameras, sensors... As a specific example, in astronomy instruments such as the Large-Synoptic Survey Telescope (LSST - [START_REF] Ivezic | [END_REF]) or the Square Kilometer Array (SKA - [Blake et al. 2004]), to only name two of them, will produce data of several terabytes and even petabytes of memory [START_REF] Longo | [END_REF]. Such unprecedented quantities require to develop scalable data processing tools being able to operate in the large-scale regime. In this work, we will focus on a specific method called Blind Source Separation (BSS), which is a key analysis tool to learn meaningful decompositions of multivalued data and which has already been successful in a wide variety of scientific fields such as audio processing [Vincent et al. 2011, Vincent et al. 2003, Ozerov & Févotte 2010, Duong et al. 2010, Févotte et al. 2009], biomedical data processing [START_REF]Removing electroencephalographic artifacts by blind source separation[END_REF], Negro et al. 2016, Poh et al. 2010] or astrophysics [Bobin et al. 2014], to only cite three of them. More specifically, the objective of this thesis is to extend a sub-domain of BSS, namely sparse BSS, to various large-scale contexts.

In the next section, we start by giving a few examples of practical BSS problems, which will be followed by a formalization of the problem. Finally, the last part of this chapter will introduce the problems faced in the large-scale setting and the structure of this work.

A.1 Four BSS examples

A.1.1 Cocktail party problem

The most well-known BSS example might be the cocktail party problem. During a party, n persons are speaking with each other. A given number m of microphones are located in the room, registering the conversations. As the guests are speaking at the same time, they create a brouhaha and the signals X 1 , X 2 , ..., X m are corresponding to different mixings of their original sentences S * 1 , S * 2 , ..., S * n . The signals captured by each microphone are different, depending on several physical effects (such as for instance, the distance of each guest to a given microphone) which are not necessarily known beforehand. As such, the mixing process leading to the mixings X 1 , X 2 , ..., X m is supposed to be unknown. The challenge of BSS is then to recover the sounds pronounced by each individual from the mixings in a blind fashion : from X 1 , X 2 , ..., X m only, we want to retrieve S * 1 , S * 2 , ..., S * n .

A.1.2 Spectroscopy : LC/MS data

The Liquid Chromatography -Mass Spectrometry (LC/MS) [START_REF] Rapin | NMF with sparse regularizations in transformed domains[END_REF]] enables to study a fluid in order to identify and quantify its constitutive chemicals. This fluid could correspond to a drink and the goal of BSS would then be to identify the spectra S * 1 , S * 2 , ..., S * n of each chemical (e.g. caffeine, sucrose, menthone...) from the LC/MS data X 1 , X 2 , ..., X m . Furthermore, if the mixing is supposed to be linear, the mixing coefficients (which are proportional to their concentrations) can also be of interest. The LC/MS data are collected within a matrix X ∈ R m×t in the following way :

-A first physical imperfect separation is performed, during which the fluid goes through a chromatography column and its chemicals are separated according to their speeds (which themselves depend on their physical properties). At each time moment 1/f e (determined beforehand by a given rate f e ), the output of the chromatography column is then sampled, giving a first imperfect time separation.

-At each time 1/f e , the samples are analysed by a mass spectrometer. Since the different chemicals have different masses, this yields a second imperfect separation in mass.

This double imperfect separation both in time and in mass gives data similar to the ones of Fig. II.1 : the time axis corresponds to the different columns of X, while the mass to charge ratio corresponds to the rows of X. It has to be emphasized that at this point, due to the imperfectness of the separations, the chemical are still mixed within the X matrix. The goal of BSS is then to enhance the separation to precisely recover the spectra of the different chemicals as rows S * 1 , S * 2 , ..., S * n of a matrix S * (cf. 3, in which it can be seen that several sources emit around 1, 5, 9 and 30), making that merely taking pictures at different wavelength values does not enable to directly separate them well. Resorting to BSS enables to estimate the elementary emissions S * i , i ∈ [1, 4], which in turns allows to study the physical mechanisms at stake in the supernovae remnant. Furthermore, BSS methods can also estimate the corresponding spectra

A * i , i ∈ [1, 4]. -SKA
The SKA is a continental-size radio-telescope that should be built in Australia and South America. It has several objectives, such as testing the general relativity, studying the large-scale structure of the cosmos, the epoch of reionization... The SKA will be the world's largest telescope, eventually comprehending thousands of dishes and up to a million low-frequency antennas2 . However, due to the wide range of frequencies in which it will operate and its size (approximately one square kilometer of collecting area), the quantities of data will be tremendous : it could indeed produce more data per day than the entire content of internet in 20113 . Therefore, while the faced BSS problems are similar in principle to the ones of the previous Chandra example, the datasets are much larger, calling for corresponding adequate BSS methods that will enable to draw as much information as possible from the data. As an example, datasets with up to m = 10000 wavelength bands and up to t = 10 9 pixels will potentially have to be considered.

A.1.4 Show-through removal

The show-through effect appears during the scanning process of documents. The issue at hand is that the backside of the document is visible on the front side, for instance due to a too high transparency of the paper [Merrikh-Bayat et al. 2011]. Thus, when reading the content of one side of the paper, the reader also sees the writing on the back. This can be written has a BSS problem with two observations X 1 and X 2 , each corresponding to the scan of one side of the paper on which the show-through effect is visible, and two sources S * 1 and S * 2 , each corresponding to an image of what is actually written on each side of the paper. The goal of BSS is to retrieve the sources, that is to remove the contribution of the other side of the paper that deteriorates the reading. See Figure II.4 for a concrete example of show-through.

A.2 Mixing model

It is now time to introduce slightly more mathematically the BSS problem, by focusing on two different possible kinds of mixings : the linear and the non-linear models.

(a) X 1 (b) X 5 (c) X 9 (d) X 30 (e) Synchrotron (f) Thermal (g) Fe 1 (h) Fe 2
Wavelength value (arbitrary units) Energy (arbitrary units) (i) A *

A.2.1 Linear model

The linear model is largely the most widespread one and has led to a large variety of works [Comon & Jutten 2010]. In these, the m observations X 1 , X 2 , ..., X m are supposed to be the linear combinations of n sources S * 1 , S * 2 , ..., S * n , each of them having t samples :

∀i ∈ [1, m], X i = n j=1 a * ij S * j (II.1)
The linear BSS model can be written in matrix form : the observations X i ∈ R t (with i ∈ [1, m]) are grouped as rows of a matrix X (of size m × t) and similarly the sources are stacked as a matrix S * (of size n × t). The model then becomes :

X = A * S * + N (II.2)
with A * (of size m × n) the mixing matrix containing the linear mixing coefficients a * ij . The matrix N (of size m × t) enables to take into account some slight deviations from the ideal linear model, such as the ones induced for instance by some noise in the measurements. Thus, N is called the noise matrix. In the current work, we will further focus exclusively on the over-determined case 4 in which n ≤ m. The objective of linear source separation is to retrieve from the sole knowledge of X the sources S * (up to limited indeterminacies, cf. below), as well as in the present blind setting the linear mixing coefficients A * . As such, BSS is a matrix factorization problem in which we aim to recover matrices having physical meanings, in contrast for instance to dictionary learning [START_REF] Mairal | [END_REF]. Written in this form, BSS is however an ill-posed problem as we are looking for a small subset of physical solutions among the infinity of possible ones. Therefore, several families of algorithms have emerged, introducing some regularization. Each of them imposes a different kind of prior knowledge on the sources, aiming to help to discriminate the desired solutions. One can cite the Independent Component Analysis (ICA - [Comon & Jutten 2010]), the Non-negative Matrix Factorization (NMF - [START_REF] Gillis | [END_REF]) and the Sparse Matrix Factorization (SMF -[Zibulevsky 2003]) families. Each one will be detailed in section A.3. A natural question is then the identifiability issue : using the previous priors on the sources, can we really hope to recover them ? Starting with the ICA methods, such questions have been well studied in the past and it has been for instance shown that the independence of the sources -in the absence of noise and Gaussian sources -enables to recover them up to a mere scaling and permutation indeterminacy (Darmois theorem [Darmois 1953]). Concerning the sparsity prior, conditions for recovery up to the same indeterminacies are studied in [START_REF] Gribonval | Dictionary Identification-Sparse Matrix-Factorization via 1 -Minimization[END_REF][START_REF] Gribonval | Sparse and spurious : dictionary learning with noise and outliers[END_REF]]. In the last work, the authors furthermore studied the non-asymptotical noisy setting in the presence of outliers, with potentially over-complete dictionaries (i.e. mixing matrices in our setting). In brief, they show that A * can be recovered using the cost function II.8, which is shown to have a local minimum around A * with high probability. Among others, some of the required hypotheses are that S * must be sparse enough (and follow other assumptions), A * must be sufficiently incoherent (depending on S * sparsity level), the noise level as well as the outlier energy are limited, and the number of non-outlier samples large enough (there are also assumptions on the regularization parameters). Note that these conditions might however be slightly restrictive in realistic experiments, and might not thus be always respected during this work.

A.2.2 Non-linear model

While convenient for many problems, the linear mixing model is only an approximation which might not always hold : it is not anymore valid when using sensors with saturations or non-linearities (for instance gas [START_REF] Madrolle | [END_REF] or chemical [START_REF] Jimenez | [END_REF], Duarte & Jutten 2014] sensors), or in some specific applications (show-through removal [Merrikh-Bayat et al. 2011], hyperspectral imaging [Dobigeon et al. 2014]). It can therefore be relevant to change the BSS model to a nonlinear one :

X = f * (S * ) + N (II.3)
Where f * is an unknown non-linear function from R n×t to R m×t (where again here n ≤ m). In this work, we will consider general functions f * , by mostly (cf. Sec. E) assuming that f * is invertible and symmetrical around the origin, as well as regular enough. More detailed hypotheses, both on the mixing and the sources, are discussed in Chapter. VII. At this point, it is important to mention that non-linear BSS is much more difficult than its linear counterpart and that it might not be possible to find both f * and S * up to a simple permutation and scaling indeterminacy. Therefore, and in contrast to usual linear BSS, the separation of the sources must be distinguished from their reconstruction :

-The sources S are said to be well separated if they are estimated up to a permutation and an unknown non-linear function h that does not remix them ;

-The sources S are said to be well reconstructed if they are estimated up to a permutation and an unknown scaling factor (thus, h is here specifically a scaling function).

It is interesting to note that in ICA, and in contrast to the linear case, the independence of the sources is not anymore sufficient to guarantee their separation. In the case of sparse sources, [Ehsandoust et al. 2016] claimed the possibility to separate the sources if only one source is active for each sample. Nevertheless, sparsity does not guarantee the reconstruction. Thus, the goal of sparse non-linear BSS is only to separate the sources by estimating the underlying non-linearities, but the source reconstruction is generally not straightforward as the problem is too ill-posed.

A.3 ICA, NMF and SMF

As said above, BSS requires tackling an ill-posed unsupervised matrix factorization problem. We here detail three additional priors (corresponding ICA, NMF, SMF) that have been introduced for enabling the identification of the mixture parameters. The discussion will be restricted to the most classical linear mixing model, the non-linear one being specifically addressed in Chapter VII. All the remainder of this work will then focus on the sparsity of the sources, that already lead to enhanced separation quality in various matrix factorization problems [Bobin et al. 2008[START_REF] Zibulevsky | [END_REF], Li et al. 2006, Le Roux et al. 2015].

A.3.1 ICA

The Independent Component Analysis is probably the family of BSS algorithms for which the most extensive literature exists. It is a statistical approach that assumes the sources to be independently identically distributed and that there is at most one non-Gaussian source. Under these assumptions, it can be shown that the sources can be recovered up to a mere permutation and scaling indeterminacy [START_REF] Darmois | Analyse générale des liaisons stochastiques : étude particulière de l'analyse factorielle linéaire[END_REF],Comon 1994] provided that the mixing matrix is invertible. The main idea is that mixing the independent sources will lead towards non independent signals, since the mixtures share the same source signals. An extensive review can be found in [Comon & Jutten 2010]. Typical ICA algorithms uses whitening and dimensionality reduction as preprocessing in order to simplify the problem. This can be for instance achieved using Principal Component Analysis (PCA - [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF], Eckart & Young 1936, Jolliffe 1986]), which enables to come back to the case m = n and to restrict the search of the mixing matrix to the group of unitary matrices -while potentially reducing the noise impact. Then, the problem boils down to finding a demixing (rotation) matrix W such that the estimated sources Ŝ = WX are independent. However, enforcing directly independence is not a trivial issue. As such, several proxys have emerged. Among them, one can cite independence measures based :

-On the mutual information (using second characteristic function or the Kullback-Leibler divergence) : this is very general, and the mutual information is zero if and only if there is independence. The drawback is however that this criterion requires the joint density of the sources, which are unknown in practice and must be replaced by an estimate, which can be relatively costly in practice.

-On contrast functions : such a function is an optimization criterion such that its global maxima correspond to a separation of all the sources. However, one then need to check that all the cumulant are zero, while in practice only a few number of them are zero.

Despite its sound mathematical foundations, the ICA family suffers from several drawbacks. First, it (implicitly) requires the knowledge of the probability density functions of the sources. While approximations can be used, these can have an impact on the final solution. Second, the mutual independence of the sources might not be a valid hypothesis in several applications (cf. Planck data [Bobin et al. 2014], or hyperspectral data : as the sum of abundance fractions associated to each pixel is constant, the sources cannot be independent [Nascimento & Dias 2005]). Finally, its performances decrease much in the presence of noise, since most the approaches are assuming noiseless mixings (and since the addition of noise create some identifiability issues of the sources [Davies 2004]).

A.3.2 NMF

The second family of algorithms builds on the special case where it is known that A * and S * have only non-negative coefficients [Paatero & Tapper 1994, Lee & Seung 1999], which is often the case on real world data [START_REF] Gillis | [END_REF]. For instance, the LC/MS data and the Chandra data of Section A.1 both fulfill this condition. The non-negativity condition states that there are vectors S * j , j ∈ [1, n] such that all the data points in X can be written as non-negative linear combinations of the S * j . In contrast to ICA, it makes that no subtraction can occur, which has led to the intuitive notion that NMF finds a data representation by patch, combining parts to form a whole [START_REF] Lee | [END_REF]. Furthermore, as the S * j are non-negative, a geometric interpretation of NMF is that it looks for a simplicial cone 5 in the positive orthant containing all the data points. However, if the data values are strictly positive, there is many such simplicial cones, raising the identifiability issue. In [Donoho & Stodden 2004], the authors however show that uniqueness of the simplicial cone can hold even if the data do not fill out the positive orthant, but

The simplicial cone generated by vectors

Sj, j ∈ [1, n] is defined as Γ S = {x ∈ R t |x = n i=1 αiSi, 0 ≤ αi}
rather a proper subset of the positive orthant. In particular, they show that the model is identifiable under the pure pixel assumption, stating that each source is active alone at least once. As this condition might not always holds, NMF has also been combined with sparsity [Hoyer 2004[START_REF] Kim | [END_REF], Rapin 2014].

A.3.3 SMF

In this thesis, we will specifically focus on sparse BSS, in which the sources are assumed to be sparse (roughly speaking, they are supposed to have a large number of zero coefficients, cf. Section B.1). This family has attracted much interest during the last two decades [START_REF] Zibulevsky | [END_REF], Bronstein et al. 2005, Li et al. 2006], which has mainly been motivated by the success of sparse signal modeling for solving very large classes of inverse problems [Starck et al. 2010]. Sparse BSS will be more extensively presented in Section B. In brief, it has lead to enhanced separation quality, in particular with the Generalized Morphological Component Analysis (GMCA - [Bobin et al. 2007]) algorithm. In the framework of ICA, Efficient FastICA (EFICA) [Koldovsky et al. 2006] is a FastICA-based algorithm that is especially adapted to retrieve sources with generalized Gaussian distributions, which includes sparse sources. In the seminal paper [Zibulevsky 2003], the author also proposed a Newton-like method for ICA called Relative Newton Algorithm (RNA), which uses quasi-maximum likelihood estimation to estimate sparse sources.

In the next section, we develop the concept of sparsity and sparse BSS.

B Sparsity and sparse BSS

B.1 Sparsity

Over the last twenty years, sparse modelling has had an increased impact in various signal processing areas, such as denoising [Elad & Aharon 2006] (or more generally signal restoration), feature extraction [START_REF] Hyvarinen | [END_REF]], compression [Le Pennec & Mallat 2000] and, in our case, source separation [START_REF] Zibulevsky | [END_REF], Bobin et al. 2007]. This interest has been mainly driven by the compressed sensing theory [Candes & Tao 2004[START_REF] Donoho | [END_REF]], which can be seen as an alternative to the Shannon sampling theory (in which the signals are assumed to be frequency band-limited and not sparse). In brief, compressed sensing theory gives some support to methods assuming the sparsity to make better-posed ill-posed inverse problems. Indeed, it can be shown that (under conditions) methods that look for a sparse solution might find the exact solution to the problem at hand. The goal of this subsection is thus to briefly present the concept of sparsity (while it is of course widely used in this work, we do not aim at presenting in depth theoretical results). The presentation will mainly follow the one of [Starck et al. 2010,Mairal et al. 2014], but we also recommend [Mallat 1999,Elad 2010], while [Rapin 2014,Che- Right : approximately sparse signal : while most of the coefficients are non-zeros, many of them are close to 0, making that the signal can be well approximated by its highest amplitude coefficients. Thus, if all the coefficients between the two red lines are set to 0, the remaining samples comprehend most of the signal energy, which makes that the approximately sparse signal is decently approximated by the exactly sparse signal of the left figure .   not 2017] give a more specific overview of sparsity for BSS.

B.1.1 Definition of sparsity and sparsity in the direct domain

Generally speaking, sparsity amounts to represent a signal with as few variables as possible. More specifically, a signal s ∈ R 1×t is said to be exactly sparse if only few of its coefficients, let say k t, are non-zeros : s 0 = k. An example of such a signal is displayed in Figure II.5a. Nevertheless, in real-life the exact sparsity is a too restrictive assumption and a more interesting case is the one of approximately sparse signals (also called weakly sparse or compressible signals). In these, while s 0 t, only a small number k of samples have a high amplitude, making that only keeping them already enables a good approximation of the signal. In this sense the signal can be well approximated as ksparse. This is for instance the case if the sorted magnitude of the samples s i , i ∈ [1, t] decays according to a power low. For an example of such an approximately sparse signal, see .6 for a concrete example). Said differently, using a transform Φ which captures well the morphology of the signal s enables to concentrate the energy of s within a few number of coefficients. Many such domains have been proposed in the signal processing community, starting from the broadly known Fourier basis [Bracewell & Bracewell 1986]. However, the Fourier basis is not localized in space, making its use restricted to stationary signals. To bypass this issue, wavelet transforms have been introduced. Compared to the Fourier basis, the main advantage is that wavelets are localized both in space and frequency domains. A wavelet basis is a set of functions Φ 1 , Φ 2 ..., Φ n that are essentially dilated and shifted version of each other. Generally speaking, wavelets enable to sparsify signals that are polynomial by parts. Starting from the ones of Haar [Haar 1910], various wavelets have been designed that correspond to different signal or image geometric contents. Among them, one can cite 6 : Daubechies wavelets [Daubechies 1988], the curvelets [Candes & Donoho 2002], the contourlets [Do & Vetterli 2005], the bandlets [Le Pennec & Mallat 2005], the starlet [Starck et al. 2010], the shearlets [Easley et al. 2008]... It is important to emphasize that among all these proposed transforms, a distinction must be made between i) orthogonal (in which Φ T Φ = ΦΦ T = Id) versus other bases of wavelets ; ii) redundant (T > t) versus nonredundant transforms. While redundant transforms enables interesting properties, they are more difficult to handle since the decomposition in the transformed domain is not unique and Φ T Φ = Id -see [Rapin 2014] for a detailed account of the use of redundant representations in the case of sparse BSS. Once a transformed domain Φ has been chosen, we still have to find a sparse representation of s in this domain. For redundant transform, such a representation is not unique and it is then possible to find among all the possible representations the one that will follow a given criterion, namely here sparsity. More specifically, for a given s ∈ R 1×t the goal is to learn a representation s s Φ Φ T such that s Φ ∈ R 1×T is sparse. An option is to learn the representation s Φ from the minimization of a cost function using the synthesis formulation :

argmin s Φ ∈R 1×T h(s Φ Φ T ) + g(s Φ ) (II.4)
with a data fidelity term h(s Φ Φ T ) ∈ R + expressing the discrepency of the model with regards to the observations s, and a regularization term g(s Φ ) ∈ R + , which is used to enforce the sparsity of the representation. From this quite general model, several specific ones can be derived, depending on the choice of h and g. For instance, a natural choice of g is to choose it using the 0 norm, to count the number of nonzeros elements : g(s Φ ) = λ s Φ 0 . However, in practice using the 1 norm instead of the 0 one is common, since the 1 norm is the closest convex norm to the 0 one [Chen et al. 2001]. Concerning h, generally a squared Euclidian 2 distance is used, 6. As a side remark, instead of imposing beforehand the transform, it is further possible to learn it from the data through dictionary learning [Olshausen & Field 1996, Mairal et al. 2014, Mensch et al. 2018], making it much more adapted to the data at hand. This is however not the main focus of this thesis.

which stems from the assumption of a Gaussian noise : h(s Φ Φ T ) = 1 2 ss Φ Φ T 2 2 . These choices lead to the well known basis pursuit formulation :

argmin s Φ ∈R 1×T 1 2 s -s Φ Φ T 2 2 + λ s Φ 1 (II.5)
For other usual formulations, we refer the reader to [START_REF] Mairal | [END_REF]. We would like to recall -while this is not the main focus of this thesis -that another formulation than the synthesis one of II.4 is possible, namely the analysis formulation :

argmin s∈R 1×t h(s) + g(sΦ) (II.6)
While for orthonormal transforms Φ, synthesis and analysis formulations are equivalent, this is not anymore the case for redundant transforms [Rapin 2014].

B.2 Sparse BSS as an optimization problem

The goal of this section is to give a quick overview of the optimization problem we will have to tackle to solve the BSS problem. More details will be given in the next chapter.

B.3 General problem

Generally speaking, in this work we will aim at performing BSS through the minimization of a non-convex cost function of the form :

argmin A∈R m×n ,S∈R n×t 1 2 h(A, S) + J (A) + G(S) (II.7)
Where :

h is a data fidelity term measuring the discrepancy between the data and the mixture model. In particular, h will be assumed to be differentiable with respect to A and S, to have Lipschitz gradients and to be block multi-convex [START_REF] Xu | [END_REF], cf. Chapter III-C.

-The penalizations J and G enforce some desired properties on A and S. The conditions on such functions will be developed in Chapter III-C. As a quick example, J and G can be used in NMF to impose the non-negativity of the sources and the mixing matrix respectively. In this situation, J (A) = ι . 0 (A) and G(S) = ι . 0 (S), where ι U is the indicator function of the set U : here, it is equal to infinity if at least one coefficient of the corresponding matrix is negative, and zero otherwise.

The minimization schemes used for Eq. II.7, and in particular the use of the multi-convex structure of h will be detailed in Chapter A.

B.4 Concrete example : sparsity prior

To make things more concrete, we now derive a more specific version of Eq. II.7 in the context of sparse BSS, in which case it can be written using the following form :

argmin 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) (II.8) -The h(A, S) = 1 2 X -AS 2
F term promotes a faithful reconstruction of the data under a white Gaussian noise N assumption inducing the use of the Frobenius norm . F (see e.g. [Hamzaoui & Bobin 2018, Bobin et al. 2019] for the case of Poisson noise).

-The G(S) = R S (SΦ T S ) 1 term promotes a 1 sparsity of the sources, with denoting the Hadamard product. This work will focus on the 1 -norm, which has been shown to help the separation of the sources in the context of sparse matrix factorization [Rapin et al. 2013]. The Φ S matrix (of size T × t, with T t) corresponds to the transform enforcing the sparsity. In this thesis, it will be taken equal to either the identity (in which case the sparsity is enforced in the direct domain) or the starlet transform [START_REF] Starck | [END_REF]]. The regularization parameter R S matrix (size n × T ) controls the trade-off between the data fidelity and the sparsity terms. It can be decomposed into R S = Λ S G where Λ S (n × n) is a diagonal matrix of the regularization parameters λ 1 , λ 2 , ..., λ n and G (n×T ) is a matrix used to introduce individual penalization coefficients in the context of reweighted 1 [START_REF] Candes | [END_REF] (when no reweighting is used, G = 1 n×T is a matrix with all coefficients equal to 1).

-To avoid degenerated A and S matrices where A F → ∞ and S F → 0 due to the sparsity penalty, the last term enforces the oblique constraint, ensuring that for all i, the i th column A i of A is onto the 2 hypersphere :

J (A) = ι {∀i∈[1,n]; A i 2 2 =1} (A)
, with ι being the characteristic function of the corresponding convex set.

B.5 Sparse BSS : hypotheses made throughout this work

In the following and without other indication, we will assume that :

-Both the matrices A * and S * are full rank : this is a common assumption, especially as in ICA the matrix A * is supposed to be orthogonal ;

-The number of sources n is known.

-We are in the over-determined setting with n ≤ m : the under-determined setting could however probably be tackled using the method of Chapter V, but this is out of the scope of this work. Furthermore, m t : high values of t are required for identifiability issues [START_REF] Gribonval | Dictionary Identification-Sparse Matrix-Factorization via 1 -Minimization[END_REF].

Chapitre II. Sparse Modelling and Large-Scale BSS -The matrix A * has unitary columns : more specifically, each column lies on the 2 unit sphere, which enables to leverage the scale indeterminacy issue (cf. Sec. A.3.3).

-The sources are sparse in a transform domain Φ S and respect the morphological diversity assumption : see Chapter A-C.4.2 for a more detailed explanation, and [START_REF] Gribonval | Dictionary Identification-Sparse Matrix-Factorization via 1 -Minimization[END_REF] for identifiability in sparse matrix factorization.

C Large-scale sparse BSS and organization of the manuscript

This thesis is subdivided into four sub-problems taking their roots around the large-scale sparse BSS issue. We here detail them and shortly present some of the main results to give a quick overview of this work. The methods, results and interpretations will be more thoroughly presented in the following Chapters, after an introduction in Chapter III of the background required for understanding the sparse BSS optimization framework.

C.1 Preliminary question : avoiding relaunchs in sparse BSS and increasing robustness

This first problem is to be linked the difficulty of finding a (local or global) minimum of the cost function II.7 corresponding to a physical factorization of the data X 7 . To perform such a task, it is possible to tune the regularization parameters R S to try to highlight the most interesting critical points. Another way is to find a good initialization of the algorithm, that is if possible an initialization within the basin of attraction of a minimum corresponding to good estimate  and Ŝ.

To perform such a task, it is usual to try several initializations and R S values, and to retain the ones corresponding to the best estimates. This grid-search approach can however become costly in the large-scale context we aim at tackling in this thesis, as the computational cost of re-launching the algorithms can become prohibitive 8 . As such, methods that are robust to the initialization and enabling an automatic regularization parameter choice without any relaunch are of uttermost importance for working on practical datasets. Therefore, in Chapter IV we will aim at trying to generalize the automatic regularization parameter choice of GMCA [Bobin et al. 2007] to the more recent and mathematically grounded Proximal Alternating Linearized 7. Thus, we highlight that we are not interested by finding any critical point of Eq. II.7, but rather a specific one corresponding to the physical true underlying factor A * and S * . Note also that there are several global minima, which can however potentially correspond to estimations  and Ŝ of unequal qualities (see e.g. [Neyshabur 2017] in the context of over-parametrized machine learning). 8. Note that this study is nevertheless not fully restricted to the large-scale setting, since in some experiments practitioners might have no clue concerning the true A * and S * , making it hard to assess the quality of the estimates  and Ŝ obtained for given parameters and initializations.

Minimization (PALM) algorithm [START_REF] Bolte | [END_REF]. The proposed strategy is empirically shown to be robust to the initialization. More precisely, Chapter IV is subdivided in two main parts :

-We first conduct an empirical study showing that, for sparse BSS, a minimization of Eq. II.8 using PALM jointly with a grid-search approach for R S choice is a strategy suffering from a low efficiency and versatility. In brief, this highlights that performing sparse BSS using PALM requires a very careful and sensitive tuning of R S . Similar conclusions can be (to a more limited extent) shown concerning PALM initialization : this approach also undergo a low reliability ;

-To alleviate the cumbersome hyper-parameter tuning described above and increase the robustness to the initialization, we then rationalize a two-step strategy, with i) a warm-up stage comprehending a GMCA, yielding in practice a descent first guess of  and Ŝ ; ii) a refinement stage using a PALM : the regularization parameters R S are chosen based the warm-up stage first guess.

The quality of the method is shown on a realistic astrophysical experiment and the limitations of the proposed two-step strategy are discussed. We refer the reader to Chapter IV.

C.2 Large number of sources n

Chapter V extends the above two-step strategy to tackle mixings with a high number of sources n. A classic example of such a problem is the one of the LC/MS data developed in Section A.1.2 : if the fluid to be analysed has a high number of constitutive elements, the corresponding number of spectra and thus of sources can reach for instance up to n = 100. No sparse BSS algorithm is currently tailored for such problems, as the issue is twofold :

-Computational issue : As the number of sources increases, so do the size of the A * and S * matrices. Since state-of-art methods such as [Bobin et al. 2007] requires several inversion of such matrices, the computational burden can become high ;

-Deteriorated performances : Beyond the computational aspects, the separation quality of most BSS algorithms tends to dramatically deteriorate in the presence of a high number of sources (cf. The approach proposed in Chapter V re-uses the two-step strategy and introduces coordinate blocks of intermediate sizes r. This is in contrast to the other stateof-art methods that either use sizes r = 1 (hierarchical of deflations methods) or full-size blocks r = n (e.g. GMCA). Beyond a huge gain in computation time, the proposed block-GMCA (bGMCA) algorithm is also empirically shown to improve the separation quality for relatively small block sizes. An interpretation of such a phenomenon, as well as extensive numerical experiments, are proposed in Chapter V.

C.3 Large datasets X

The most tangible large-scale issue might the one of large datasets X ∈ R m×t . As evoked in introduction, these might for instance become the new standard in astronomy, with space telescopes such as Euclid 9 , huge radio-interferometers such as the SKA 10 (cf. Sec. A.1.3), or the telescope LSST 11 . Such devices will produce tremendous amounts of data : in order to give a rough idea of the involved quantities, one could speak of images of resolution t = 10 9 pixels and about m = 10 4 channels. In this context, sparse BSS algorithms that yielded high quality results [Bobin et al. 2007] cannot be used, due to :

-Memory limitations : The datasets are so huge that they cannot be stored at once within memory. Therefore, one of the few current solutions is to scan submatrices of X and work on each of them in an isolated way, which might deteriorate much the separation quality ;

-Time limitations : Even not considering memory issues, the heavy computations involved might largely slow down the processing time. As such, parallelized algorithms should be preferred.

Chapter VI proposes to introduce mini-batches, that is submatrices of X, within the GMCA algorithm 12 . Note that this is in contrast to the use of intermediate-size block coordinate methods described above, as these split the factors A and S but use the whole X. However, the use of mini-batches with the projected Alternating Least-Square (pALS) scheme of GMCA raises open questions, as each of the B minibatches yields a different estimate of the full A * . Therefore, we need to aggregate these estimates at each iteration of the proposed distributed-GMCA (dGMCA). In brief, we propose to benefit from the structure implied by the oblique constraint in Eq. II.8 to enhance such an aggregation through the use of a Fréchet mean. To further robustify the process, we also propose to i) use weights taking into account the assumed quality of the estimation of each source in each mini-batch ; ii) use a more robust aggregation method, reminiscent of a median on the unit hypersphere.

Beyond the expected gain in terms of computation time and memory (as each minibatch can be sent to a different node of a cluster), a striking result is that when the sources S * are very sparse, the use of the robust Fréchet mean enables to obtain better results than with full batch methods. We propose an explanation of this phenomenon in Chapter VI.

C.4 Non-linear BSS

Non-linear BSS is a difficult problem, although such a model can be required when the usual linear one is too simplistic for correctly modelling the mixing. An 9. http ://sci.esa.int/euclid/ 10. https ://www.skatelescope.org 11. https ://www.lsst.org 12. Indeed, the two-step approach requires as a warm-up stage at least a few iterations of GMCA. As such, it still needs to be able to use GMCA on large X.

example is the show-through removal presented in Section A.1.4, for which the physical processes imply a linear quadratic model rather than a mere linear model. The StackedAMCA algorithm we propose in Chapter VII to tackle non-linear sparse BSS can be envisioned as an extension of large-scale BSS. Indeed, it estimates the underlying non-linearities through a linear-by-part model comprehending L parts. As this approach requires to solve a potentially large number L of linear subproblems, it calls for efficient linear BSS algorithms. More specifically, StackedAMCA is an iterative algorithm. Each iteration alternates between a linear step and a non-linear one :

-The linear step consists in the application on the current dataset of a sparse linear BSS algorithm robust to high amplitude non-linearities, which enables to find a new part of the linear-by-part approximation ;

-The non-linear step consists in creating a new dataset in which the contribution of all the previously found linear models is canceled. This paves the way for the next iteration, in which a new linear model will be fitted on this new dataset.

The algorithm obtains better separation results than other state-of-art methods. Interestingly enough, it can further reconstruct the sources in some settings despite increased indeterminacies in non-linear BSS (cf. Sec. A.2.2). The required hypotheses for StackedAMCA to work, as well as a subset of mixings for which it can further be expected to reconstruct the sources, are discussed in Chapter VII.

Chapitre III

Optimization frameworks for sparse BSS

A Reminder about sparse BSS cost function and outline of the chapter

As mentioned in the previous chapter, we will aim at performing BSS through the minimization of a cost function of the form :

argmin A∈R m×n ,S∈R n×t 1 2 h(A, S) + J (A) + G(S) (III.1)
Such a minimization is not easy to perform due to several reasons :

-While the term h(A, S) is assumed to be smooth with a Lipschitz gradient, calling for simple optimization methods such as gradient descent, it is not necessarily the case of the J (A) and G(S) terms (cf. e.g. Sec. B.4 explaining the case of sparse sources and the use of the 1 -norm, which is not differentiable in 0). Therefore, we will need to resort to non-smooth optimization algorithms and more specifically to proximal algorithms, which will be explained in section B.

-The problem is non-convex in A and S. Nevertheless, as mentioned in Section B.3, it is multi-convex [START_REF] Xu | [END_REF], calling for the use of specific algorithms that will be detailed in sections C.3 (algorithms aiming at truly minimizing Eq. II.7), and C.4 (heuristics enabling enhanced practical results).

We assume the reader to have a basic knowledge of optimization tools, and refer to [START_REF] Boyd | [END_REF] for generic notions.

B Proximal operators and proximal algorithms

Proximal algorithms enable to solve convex non-smooth optimization problems. They have been extensively used over the last decades in the signal processing community [START_REF] Combettes | [END_REF], Jenatton et al. 2010, Bolte et al. 2010, Jezierska et al. 2012,Bolte et al. 2014,Chouzenoux et al. 2016]. Roughly speaking, they operate at a high level of abstraction, since their base operation uses the proximal operator of a function, which computation involves itself a (potentially simple) convex optimization problem. In this section, we will first explain the notion of proximal operator and then their use within proximal algorithms. Our presentation will mainly follow the ones from [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF], Parikh et al. 2014].

B.1 Proximal operators

B.1.1 Definition

We will call a function f : R t → R ∪ {∞} proximable if it is a closed proper convex function. For such functions, the proximal operator of f is defined as :

prox f (u) = argmin y∈R t f (y) + 1 2 u -y 2 2 (III.2)
with .

2 the usual Euclidian norm for vectors. Since the minimized function is strongly convex and not infinite everywhere, it has a unique minimizer. Several interpretations of the proximal operators are interesting for their understanding :

-A natural one is to see proximal operators as a local minimizer of f : the function minimized within the operator is the sum of f and a term that penalizes high distances between the current point u and y. More specifically, if u m is a minimum of f , then prox f (u m ) = u m , which is called the fixed point property. This is one of the bases for understanding proximal algorithms and, combined to the firm non-expansiveness of proximal operators, leads to one of the most basic method, namely the proximal point algorithm (cf. Sec. B.3.1).

-Proximal operators can be interpreted as gradient steps for minimizing f or a function related to f . For instance, introducing the Moreau envelope of f :

M f (u) = inf y∈R t (f (y)+ 1 2 u -y 2
2 ), it can be shown that the proximal operator of f can be written as :

prox f (u) = u -∇M f (u) (III.3)
Therefore, prox f can be viewed as a gradient step with step 1 (generalization to other steps are possible) for minimizing M f , which has the same minimizers as f . This leads to a natural interpretation for minimization algorithms.

-Proximal operators can also be seen as a generalization of projections. This is explicit if f is the indicator function of a closed nonempty convex set U : f (y) = i .∈U (y). Then the proximal operator of f reduces to the Euclidian projection onto U :

prox f (u) = Π U (u) = argmin y∈U u -y 2 (III.4)
This interpretation enables to understand the use of some proximal operators in proximal algorithms as a projection onto some constraint set.

While the requirement of computing the proximal operators through a minimization problem might seem intricate, in practice many usual functions have explicit or easy-to-compute proximal operators. The role of the following subsection is to present some of them which will be of use in the remaining of this work. But before, we highlight a last generic property of proximal operators that will be implicitly used throughout this whole work. Indeed, we will in practice be applying proximal operators on matrices, and not on vectors as defined above. Fortunately, this is not an issue since our constraints will be separable, enabling to use the separable sum property. More specifically, if f is separable such that f (y) = K i=1 f i (y i ), K ∈ N * , then :

(prox f (u)) i = prox f i (u i ) (III.5)
Where (prox f (u)) i is the i th line of the vector prox f (u). Said differently, instead of working directly on matrices, we can come back to the vector case (e.g. with the oblique constraint on A, cf. below) or even to the scalar case (e.g. with the sparsity constraint on S).

B.2 Examples of proximal operators used in this work

This subsection gives a few examples of proximal operators, by focusing on the constraints J and G that will be used in this work.

1 -Penalizations G for the sources S and corresponding proximal operators :

-1 sparsity constraint in some transformed domain : In this case, the sparsity constraint on S is enforced with a 1 -norm penalization (which has led to enhanced separation quality in BSS [Rapin et al. 2013]) :

G(S) = R S (SΦ T S ) 1 (III.6)
where the sign denotes the Hadamard product. The matrix R S (of same size as SΦ T S ) contains individual regularization parameters for each coefficient. The Φ S is a transform, supposed orthogonal in the following, into a domain in which S can be sparsely represented. The proximal operator of G is then explicit and corresponds to the soft-thresholding operator with threshold R S -which we shall denote S R S , cf. Appendix A for the definition -applied in the transformed domain : prox G (S) = S R S (SΦ S T )Φ S -Non-negativity in the direct domain : here, all coefficients in S must be non-negative, which can be written using the characteristic function of the positive orthant

K + = {S ∈ R n×t ; ∀i ∈ [1, n], j ∈ [1, t], S j i ≥ 0} : G(S) = ι K + (S) (III.7)
Following the previous subsection B.1.1, the corresponding proximal operator is then the projection on the positive orthant Π K + , which is the identity for non-negative coefficients and 0 for the negative ones (cf. Appendix A for the exact definition).

-Non-negativity in the direct domain and 1 sparsity constraint in some transformed domain : due to the non-negativity constraint, all coefficients in S must be non-negative in the direct domain in addition to the sparsity constraint in a transformed domain Φ S . It can be formulated as :

G(S) = R S SΦ S T 1 + ι K + (S) (III.8)
The difficulty is to enforce at the same time two constraints in two different domains, and the corresponding proximal operator of G is not explicit. It can either be roughly approximated by composing the proximal operators of the individual penalizations to produce a cheap update : prox G (S) = S R S (Π K + (S)Φ S T )Φ S , or computed accurately using the Generalized Forward-Backward splitting algorithm [START_REF] Raguet | [END_REF] (cf. Sec. B.3.3).

2 -Penalizations J for the mixing matrix A :

-Oblique constraint : the columns of A are constrained to lie on the 2 hyper-sphere. The corresponding set is

O = {A ∈ R m×n ; ∀j ∈ [1, n], A j 2 2 = 1}.
This constraint is used to avoid degenerated  and Ŝ results. More specifically, J can be written as :

J (A) = ι O (A) (III.9)
Following this constraint, the proximal operator is prox J (A) = Π O (A), the projection on the 2 unit hypersphere of each column of A (cf. Appendix A).

-Non-negativity and oblique constraint : Adding the non-negativity constraint to the oblique constraint on A reads :

J (A) = ι O (A) + ι K + (A) (III.10)
The proximal operator can be shown to be the composition of the proximal operator corresponding to non-negativity Π K + followed by

Π O : prox J (A) = Π O (Π K + (A)).

B.3 Proximal algorithms

Now that the proximal operators have been introduced, we can explain how to use them for solving convex (non-smooth) optimization problems. To do that, we will detail several so-called proximal algorithms. While the first one is mostly described for the sake of clarity and introducing the concept of proximal algorithms, the two following ones will be cornerstones of this work1 .

B.3.1 Proximal point algorithm

Let f : R t → R be a closed proper convex function. The problem here is to find a minimizer of f : argmin

u∈R t f (u) (III.11)
To perform such a task, the most simple proximal algorithm is the proximal point algorithm, which iterations are described by :

û(l+1) = prox λf (û (l) ) (III.12)
where û(l) it the l th iteration of the algorithm and λ > 0. Then, if f has a minimum, the algorithm converges to it, which is linked to the fact that the proximal operator is a non-expansive operator and a local minimizer of f . However, this simple method does not have many applications, since it is restricted to the cases where f is difficult to minimize but f plus a quadratic is simple to minimize. Therefore, more advanced minimization schemes using splitting methods have been introduced.

B.3.2 Forward backward splitting

The problem at hand is here to minimize :

argmin u∈R t h(u) + J (u) (III.13)
with h : R t → R and J : R t → R ∪ {∞} closed proper convex and h differentiable.

As an example of such a problem, note that the BSS problem of Eq. II.8 when performing the minimization over only one variable (e.g. fixing A and performing the minimization over S only), is an instance (generalizing the discussion to matrices) of such a minimization : argmin

S∈R n×t 1 2 X -AS 2 F + R S (SΦ T S ) 1 (III.14)
In accordance to its name, the Forward-Backward Splitting (FBS - [START_REF] Combettes | [END_REF]) method splits the objective into two terms, using the differentiability of h. The corresponding iterations write as :

û(l+1) = prox ηJ (û (l) -η∇h(û (l) )) (III.15)
where η > 0 is a step size and ∇h is the gradient of h. The name of the algorithm thus stems from the fact that it alternates between a forward / explicit gradient step and a backward / implicit step, namely the application of the proximal operator (special cases include : i) when J = 0, the usual gradient descent ; ii) when h = 0, the previous proximal point method). Provided that h is L-Lipchitz, this methods converges when η ∈ (0, 2/L). It has also to be emphasized that while the FBS algorithm has a convergence rate of O(1/l), it can be accelerated to O(1/l 2 ) using a linear combination of the previous estimates, which is known as the FISTA algorithm [Beck & Teboulle 2009].

While the FBS already found many applications and will be extensively used throughout this thesis, it is computationaly efficient only if the proximal operator of J is easy to compute (or, even better, explicit). If it is not the case, but if J is the sum of proximable terms which proximal operators are easy to compute, the following Generalized Forward Backward Splitting (GFBS) algorithm can be used.

B.3.3 Generalized forward backward splitting

While the FBS algorithm is restricted to the sum of one differentiable and one non-differentiable term, the GFBS [START_REF] Raguet | [END_REF] generalizes to a higher number of non-differentiable terms :

argmin u∈R t h(u) + K i=1 J (i) (u) (III.16)
with K ≥ 1 and J (1) , ..., J (K) some lower semi-continuous proper and convex function from R t to R ∪ {∞}. As an example of such a problem, looking for the sources with a fixed A in BSS when a sparsity and a non-negativity prior are both enforced leads to the following problem :

argmin

S∈R n×t 1 2 X -AS 2 F + R S (SΦ T S ) 1 + ι K + (S) (III.17)
where thus K = 2 and J (1) and J (2) would correspond to the second and third terms respectively. In this example, while the proximal operator of J (1) and J (2) can be computed easily (at least, provided that the transform Φ S is orthogonal), it is not the case of the proximal operator of the sum J (1) + J (2) , which is why a GFBS algorithm is useful.

In brief, the GFBS algorithm uses at each iteration the proximal operator of each individual constraint J (1) , ..., J (K) . Each iteration can be written as :

for i = 1..K do z (i) = z (i) + µ prox η ω i J (i) (2û (l) -z (i) -η∇h(û (l) )) -û(l) û(l+1) ← b i=1 ω i z (i)
where the ω i are weights in the sum (thus, ∀i ∈ [1, K], ω i ∈ (0, 1) and K i=1 ω i = 1), η ∈ (0, 2L) is the gradient step, bounded by the Lipschitz constant L of h, and µ ∈ (0, min(1.5, 0.5 + 1 Lη )). Thus, from the sole knowledge of the proximal operator of the J (i) 's, the GFBS enables to minimize problem III.16.

C Multi-convex optimization for sparse BSS

C.1 Problem statement

The previous section has enlighten how to handle the non-smoothness of matrix factorization problems by using proximal algorithms. We now aim at reviewing some classical sparse BSS optimization methods. Such algorithms however need to tackle the remaining non-convexity issue of the cost function II.7. Generally speaking, in this part we aim at minimizing cost functions of the following form : argmin

U (1) ∈R n 1 ×t 1 ,...,U (K) ∈R n K ×t K h(U (1) , ..., U (K) ) + K i=1 J (i) (U (i) ) (III.18)
with the following assumptions :

-The function h :

R n 1 ×t 1 ×R n 2 ×t 2 ...×R n K ×t K → R is assumed to be block multi- convex, that is for all i ∈ [1, K], h is a convex function of U (i)
while all the other blocks are fixed. We will further write h i the corresponding convex function over one block :

h (i) : R n i ×t i , U → h(U (1) , ..., U (i-1) , U, U (i+1) ..., U (K) ) for any arbitrary U (1) ∈ R n 1 ×t 1 , ..., U (i-1) ∈ R n i-1 ×t i-1 , U (i+1) ∈ R n i+1 ×t i+1 , ..., U (K) ∈ R n K ×t K .
In the following, each h (i) is further assumed to be differentiable and its gradient to be L i -Lipschitz ;

-For all i ∈ [1, K], the extended-valued functions J (i) : R n i ×t i → R ∪ {∞} are assumed to be convex, proper and lower-continuous. While these functions are used to enforce only individual penalizations over the U (i) , we however point out that they can be non-smooth.

Thus, the algorithms we will present in this section encompass a wider setting that the actual BSS problem of II.7, which is a specific case of Eq. III.18 with K = 2 and U (1) = A, U (2) = S.

C.2 Overview

All the algorithms we present in this section are part of the Gauss-Seidel family. More specifically, they use the multi-convex structure of the cost function (III.18) by performing an (approximate) minimization over each block U (1) , ..., U (K) . We will furthermore distinguish two sub-families of algorithms :

-Algorithms trying to minimize exactly problem (III.18) :

The multi-convex problem described by Eq. (III.18) can be tackled in its exact form using one of the following three main algorithms : Block-Coordinate Descent (BCD - [Tseng 2001]), Proximal Alternating Linearized Minimization (PALM - [START_REF] Bolte | [END_REF]) and Proximal Block Coordinate (PBC - [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems : An approach based on the Kurdyka-Łojasiewicz inequality[END_REF]). Interestingly, these algorithms can precisely converge to a critical point of (III.18).

-Algorithms minimizing an approximation of problem (III.18) : Among them, the Projected Alternating Least Squares (pALS) was first introduced in the context of NMF [Paatero & Tapper 1994]. Contrary to the previous algorithms, it does not truly minimize (III.18), but rather an approximation of it. In the context of sparse BSS, we will further discuss one of pALS extensions which has known a wild success : the Generalized Morphological Component Analysis (GMCA [Bobin et al. 2007]). GMCA is based on the pALS framework, but uses an automatic decreasing hyper-parameter strategy along the iterations that yields much more robustness in practice. In addition, while the convergence of GMCA is not guaranteed because it is built on pALS scheme, in practice this heuristic 2 helps it to numerically stabilize.

C.3 Algorithms aiming at finding a (local) minimum

We first review the algorithms that aim at minimizing exactly problem III.18. They all have a strong mathematical background, and in particular they are proved to converge under our assumptions [START_REF] Xu | [END_REF] (cf. Appendix B for more details concerning this topic). It is however important to emphasize that due to the nonconvexity of the problem, the minimizers yielded by the different algorithms might be different.

C.3.1 BCD

The BCD algorithm [Tseng 2001] was one of the first algorithms to be able to tackle Eq. III.18. It uses the multi-convex structure of the problem by solving exactly and cyclicly convex subproblems of the form argmin

U (i) ∈R n i ×t i h (i) (U (i) ) + J (i) (U (i) ) (III.19)
Each of these subproblems can then be handled for instance by the usual proximal algorithms presented in section B.3, depending on the difficulty to compute the proximal operator of J (i) (in particular, if such a proximal operator is explicit, a FBS [START_REF] Combettes | [END_REF] or one of its accelerated version is a straightforward choice). BCD is summarized in Algorithm 1, where h

(l) (i) is a short-hand for h (l) (i) : R n i ×t i , U → h(U (l) (1) , ..., U (l) (i-1) , U, U (l-1) (i+1) ..., U (l-1) (K) ), that is h (i)
with the fixed matrices chosen as the previously updated matrices. While the convergence conditions are in practice verified in sparse BSS, more generally the minimum in each BCD step must be uniquely attained to prove the convergence [Zangwill 1969], otherwise the method may cycle indefinitely [Powell 1973]. Assuming a strict convexity of each subproblem III.19 then enables to prove that every limit point of the sequence of iterates is a critical point of III.18.

2. By heuristic, we here mean an approximate method enabling to obtain non-optimal but still decent results, while alleviating a difficulty of the initial problem (e.g. time computation issues, difficult hyper-parameter choice, need to perform relaunches with different initializations...). As such, the decreasing parameter choice of GMCA is an heuristic in that the parameters should rigorously stay fixed during the whole algorithm to minimize Eq. II.8 : we shall see in Chapter IV that the corresponding estimates  and Ŝ are good while sub-optimal compared to what could be obtained with optimal fixed regularization parameters. This heuristic however alleviate the cumbersome search of regularization parameters.

Algorithm 1 BCD function BCD( Û(0) (1) , ..., Û(0) (K) ) Û(0) (1) , ..., Û(0) (K) is the initialization while not converged do for i = 1...K do Û(l+1) (i) = argmin Û(i) h (l) (i) ( Û(i) ) + J (i) ( Û(i) ) end for l ← l + 1 end while return Û(l) (1) , ..., Û(l) (K) end function C.3.2 PBC
Removing this strict convexity assumption can be done by relaxing the minimization of each subproblem III.19 with a proximal term [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems : An approach based on the Kurdyka-Łojasiewicz inequality[END_REF], yielding new subproblems that each corresponds to a local minimization : argmin

U (i) ∈R n i ×t i h (l) (i) (U (i) ) + J (i) (U (i) ) + 1 2L (l) (i) U (i) -U (l) (i) 2 (III.20) With the L (l) (i)
positive numbers (that are finite). This is, by definition of proximal operators, equivalent to : prox

L (i) ×(h (l) (i) +J (i) ) (U (l) (i) ) (III.21)
Hence the name of the Proximal Block Coordinate -PBC -algorithm. This leads to the iterative scheme of Algorithm 2.

Algorithm 2 PBC function PBC( Û(0) (1) , ..., Û(0) (K) ) Û(0) (1) , ..., Û(0) (K) is the initialization while not converged do for i = 1...K do Û(l+1) (i) = prox L (l) (i) ×(h (l) (i) +J (i) ) ( Û(l) (i) ) end for l ← l + 1 end while return Û(l) (1) , ..., Û(l) (K)
end function

While this algorithm has interesting properties (in particular, it is proved to converge under mild conditions and might be intuitively faster than BCD [Chenot 2017]), it has not been extensively used in sparse BSS. This might be linked to the fact that it requires the proximal operator of h (l) (i) + J (i) to be explicit to yield computationally cheap updates. If this is not the case, such a proximal operator must be computed using a subroutine (e.g. a FBS), which however leads to the non-trivial issue of computational error accumulations in each step. Thus, it has been claimed in [START_REF] Bolte | [END_REF]] that such a scheme is mainly "conceptual", calling instead for the PALM algorithm. Therefore, the PBC will not be studied in this work.

C.3.3 PALM

PALM has been introduced in [START_REF] Bolte | [END_REF], Xu & Yin 2014] and has driven an intensive research with a large number of applications and extensions3 in the scope of matrix factorization [START_REF] Chouzenoux | [END_REF], Chouzenoux et al. 2016]. It can be seen as the merging of PBC with the FBS algorithm, or said differently, as an alternating minimization approach for the FBS. Compared to BCD, instead of fully minimizing the subproblems III.19, PALM minimizes a local regularization of the Gauss-Seidel scheme. However, in contrast to PBC, it uses a proximal linearization of each subproblem. Therefore, it avoids the entirely implicit step required by PBC. More specifically, the minimization performed by PALM is : argmin

U (i) ∈R n i ×t i h (l) (i) (U (l) (i) ) + ∇h (l) (i) (U (l) (i) )|U (i) -U (l) (i) + J (i) (U (i) ) + L (l) (i) 2 U (i) -U (l) (i) 2
(III.22) where L (i) is a constant that we will take equal to the Lipschitz constant of ∇h

(l) (i)
in the following. Compared to BCD, the advantage is that due to the linearization, U (i) do not appear anymore in the first term. Thus, the update can be rewritten using only the proximal operator of J (i) : prox

J (i) L (l) (i)   U (l) (i) - 1 L (l) (i) ∇h (l) (i) U (l) (i)   (III.23)
Therefore, for each block of coordinates, PALM performs one gradient step on the smooth part, while a proximal step is taken on the non-smooth part. Thus, provided that the proximal operators of the individual regularization terms J (i) are explicit, PALM iterations have a low computational cost. In particular, it can prove in some settings to be faster than BCD, since it does not require to perform a full minimization of each subproblem at all iterations [START_REF] Xu | [END_REF]. The whole scheme is summarized in Algorithm 3.

Algorithm 3 PALM function PALM( Û(0) (1) , ..., Û(0) (K) ) Û(0) (1) , ..., Û(0) (K) is the initialization while not converged do for i = 1...K do Û(l+1) (i) = prox J (i) L (i) Û(l) (i) -1 L (i) ∇h (l) (i) Û(l) (i) end for l ← l + 1 end while return Û(l) (1) , ..., Û(l) (K)
end function

Example of application of PALM : sparse BSS Since PALM will be extensively used in this work, we here derive it for the specific case of sparse BSS corresponding to Eq. (II.8), which is done using the examples of proximal operators given in section B.2 : l) has not reached the desired value, iterate over (l) :

PALM( Â(0) , Ŝ(0) ) Requires : X, R S While the stopping criterion ∆ (
1 -Update of S using the current version of Â(l-1) :

S = Ŝ(l-1) - γ Â(l-1)T Â(l-1) 2 Â(l-1)T ( Â(l-1) Ŝ(l-1) -X) (III.24) Ŝ(l) = S γR S Â(l-1)T Â(l-1) 2 ( SΦ T S )Φ S (III.25)
2 -Update of A using the current version of Ŝ(l) :

à = Â(l-1) - δ Ŝ(l) Ŝ(l)T 2 ( Â(l-1) Ŝ(l) -X) Ŝ(l)T (III.26) Â(l) = Π . 2 =1 ( Ã) (III.27) 3 -Update stopping criterion ∆ (l) = min j∈[1,n] Â(l)j Â(l)j F , Â(l-1)j Â(l-1)j F
The adaptation of PALM to our cost function is detailed below :

-Update of S : the term Â(l-1)T ( Â(l-1) Ŝ(l-1) -X) is the gradient of the data fidelity term with respect to S and a Lipschitz modulus upper bound can be chosen as Â(l-1)T Â(l-1)

2

, where . 2 is the spectral norm 4 . The parameter 4. More specifically, if U is a matrix, x a vector and .

2 is the 2 norm for vectors, the . 2 γ is chosen in the range (0, 1), which is required to ensure the convergence 5 .

-Update of A : the term ( Â(l-1) Ŝ(l) -X) Ŝ(l)T is the gradient of the data fidelity term with respect to S and a Lipschitz modulus upper bound can be chosen as Ŝ(l) Ŝ(l)T

2

. The parameter δ is chosen in the range (0, 1), which is required to ensure the convergence.

-Stopping criterion : the stopping criterion ∆ (l) was chosen in this thesis as the cosine of the maximum angle between the columns of Â(l) and that of the previous estimate of the mixing matrix Â(l-1) . The algorithm stops when ∆ (l) becomes higher than a threshold τ fixed by the user, that is when the changes in A become very small.

N.B. : while the blocks here were chosen as U (1) = A and U (2) = S, other choices of coordinates are possible, which will be detailed in chapter V.

C.4 Algorithms minimizing an approximation of the cost function

We will now review some approximate algorithms. While such algorithms can both be cheaper and easier to use than exact ones and obtain good practical results, they however suffer from much less mathematical support.

C.4.1 pALS

The projected Alternating Least Square has first been introduced in the context of NMF [START_REF] Berry | [END_REF]] by [Paatero & Tapper 1994]. 6 As all the previous algorithms, it performs a cyclic update over all the coordinates. However, it does not perform a true minimization of each subproblem III.19 as in BCD, but rather uses a rough approximation : for a given i ∈ [1, K], the subproblem is first minimized without taking into account the constraint J (i) ; then the solution is projected through the proximal operator of J (i) . The whole procedure is described in Algorithm 4. While this algorithm can be interesting in many practical setting, especially when the minimization of h

(l) (i) (U (i)
) is explicit as well as the proximal operator of J (i) , it suffers from a lack of mathematical grounding. In particular, it is neither proved to converge nor to minimize the cost function III.18 in general. Indeed, even in the historical case of NMF, this algorithm can increase the cost function when performing the minimization over one fixed variable, which precludes any convergence guarantees in the BCD framework [Kim et al. 2008

] (cf. Appendix B).

induced matrix norm is defined as :

U 2 = sup x =0 Ux 2 x 2 (III.28)
5. Note that such a parameter might influence the quality of the solution found by the algorithm, due to the non-convexity of the cost function. Such a study is however beyond the scope of this work.

6. While it historically started with the standard projected least-square approach, h being a quadratic loss, we will present it here in a more general setting for the sake of continuity with the previous subsection.

Algorithm 4 pALS

function pALS( Û(0) (1) , ..., Û(0) (K) ) Û(0) (1) , ..., Û(0) (K)
is the initialization while not stabilized or maximum number of iterations not reached do for i = 1...K do Ũ(l+1)

(i) = argmin Û(i) h (l) (i) ( Û(i) ) Û(l+1) (i) = prox J (i) Ũ(l+1) (i) end for l ← l + 1 end while return Û(l) (1) , ..., Û(l) (K)
end function

Despite this lack of mathematical guarantees, its simplicity has however made pALS successful in the NMF community. Indeed, when the constraints J (i) are fairly simple, it can give good practical results. Concerning sparse BSS, pALS has known a huge success in the context of the GMCA algorithm. However, when the constraints become more complicated than mere NMF or sparsity in the direct domain (such as for instance trying to combine both non negativity in the direct domain and sparsity in another domain using redundant transforms), its separation performances can deteriorate in comparison to a true BCD [Rapin et al. 2013].

pALS in sparse BSS Since the pALS is the basis of the GMCA algorithm [Bobin et al. 2007], which is the core of this PhD (with PALM), we here detail it in the context of sparse BSS. However and as explained above, contrary to the PALM algorithm of section C.3.3, pALS does not truly look for a minimizer of II.8 but rather for a minimizer of an approximation of the cost function. In that, the sparsity parameters of pALS do not fully correspond to R S , the ones of II.8, and we shall denote them as M S . The algorithm is then given by (for the sake of simplicity, Φ S is supposed to be the identity matrix) :

pALS( Â(0) , Ŝ(0) )
Requires X, M S While not stabilized or maximum number of iterations not reached, iterate over (l) :

1 -S is updated using the current Â(l-1) .

Ŝ(l) = S M S Â(l-1) † X (III.29) 2 -A is updated using the current Ŝ(l) : Â(l) = Π . 2 =1 X Ŝ(l) † (III.30)
In this algorithm, the minimization of the quadratic data fidelity term of II.8 over one of the two matrices A, S is explicit and is performed using the Moore-Penrose pseudo-inverse [Ben-Israel & Greville 2003] of the fixed matrix, denoted as † . As a side remark, for this specific application of pALS to sparse BSS and in the case where at each iteration l, Â(l) and Ŝ(l) are orthogonal, the pseudo-inverse of the matrices is equal to the transpose. Thus, the PALM algorithm trivially reduces to the pALS, and as such in this case the pALS both truly minimizes the cost function II.8 and is proved to converge (and

M S = R S ). 7
C.4.2 GMCA as an optimized pALS for sparse BSS GMCA is not a new optimization framework, but rather an enhancement of the pALS algorithm in the sparse BSS case of II.7. As explained in the previous subsection, while the pALS framework is appealing due to its simplicity and interpretability, it nevertheless suffers from several issues (both theoretical and practical, such as a lack of robustness when the thresholds are too low [Chenot 2017]). The GMCA algorithm [Bobin et al. 2007] enables to both assuage such issues and to yield good practical separation results [Bobin et al. 2008, Bobin et al. 2015], while alleviating the cumbersome hyperparameter choice. More specifically, while GMCA is very similar to the algorithm of section C.4.1, its main strength is to propose an automatic adaptive parameter choice for M S , which is based on a fixed point argument and further enables to benefit from the morphological diversity assumption [Bobin et al. 2007]. In brief (the required qualities needed for a good hyperparameter choice will be more extensively discussed in chapter IV), such an automatic choice enables pALS to be more robust to the noise, more reliable (insensitive to the initialization) and more accurate. In the following, we will first state the morphological diversity principle and then the GMCA automatic hyper-parameter choice.

Morphological diversity

The morphological diversity dates back to mono-channel component separation, and has first been proposed in the context of morphological component analysis (MCA) [Starck et al. 2010]. In brief, in MCA the mixing x ∈ R t is assumed to be a linear combination of morphological components

x * i ∈ R t : x = n i=1 x * i .
Each of this components x * i is assumed to be sparsest in its own dictionary Φ i , all the Φ i being different, which corresponds to the assumption that all the x * i have different morphologies (said differently, different kinds of geometrical features). For instance, in the case n = 2, x * 1 could be broadly distributed while x * 2 could resemble punctual sources : thus, x * 1 and x * 2 would be sparsest in different dictionaries, encoding different kind of geometrical features. Using such an information 7. Note that instead of choosing the blocks as A and S but rather splitting the problem into 2n blocks corresponding to the columns of A and the rows of S, the update is explicit and pALS also becomes exact (since ∀i ∈ [0, n],

A i F = 1 and A i † = A i T ).
then enables to recover x * 1 , ..., x * n from x by looking for a minimizer of : argmin

x 1 ∈R t ,...,xn∈R t x - n i=1 x i 2 2 + n i=1 λ i x i Φ T i 0 (III.31)
On the other hand, this principle has been extended for sparse BSS to the case of components having the same kind of geometrical features, that is, that are sparse in a same dictionary Φ S . In that case, sparsity can still be used to discriminate the components, which are now the sources S * i . As the S * i are sparse in Φ S , the information is encoded in a small number of significant entries. Since the S * i are however different from each other, these significant entries are also likely to differ through the position at which they are active. Therefore, the discrimination between the sources can rely on their most significant coefficients in Φ S (cf. e.g. It is important to point out that many practical sources satisfies the morphological diversity hypothesis. In particular, i.i.d. sources drawn according to a Bernoulli-Gaussian, a Laplacian, a Generalized Gaussian with parameter α ≤ 1 or a strongly sub-gaussian distribution respect such an assumption. In the remaining of this subsection, to explain the GMCA algorithm, we assume for the sake of clarity that Φ T S = I d , that is the sources are sparse in the direct domain. The generalization to sparsity in transformed domain is relatively straightforward. 2) having its highest amplitude samples lying on the axes. This is typical of sources respecting the morphological diversity assumption, since the highest amplitude samples supports are (almost) disjoint. On the other hand, any mixing by a non-trivial A * will break such a shape since the coefficients of the largest amplitude samples will be simultaneously large (cf. Fig. III.2). Such a loss of compressibility can be measured through the 1 norm : the samples of the sources are enclosed into a 1 norm ball of smaller radius than the ones of the mixing. Thus, a good unmixing matrix  should ensure that the unmixing S =  † X =  † A * X lies into a 1 norm of smallest radius.

Another point of view on sparse BSS : a geometrical interpretation

GMCA : a pALS scheme enhanced with an heuristic threshold choice We can now explain how the GMCA algorithm takes advantage from the morphological diversity to perform sparse BSS. We will then detail the algorithm through a second interpretation in terms of a fixed point condition and noise removal.

A decreasing threshold choice enabling to explicitly benefit from the morphological diversity assumption As stated just above, the sparse BSS 1 and S * 2 . Right, down : Scatter plot of the wavelet coefficients of the two mixings. The green square corresponds to the 1 ball of minimum radius enclosing all the mixing wavelet coefficients. As expected, its radius is larger than the one of the 1 ball found with the source wavelet coefficients, in red.

S * 1 X 1 S * 2 X 2 (S * Φ T S ) 1 (S * Φ T S ) 2 (XΦ T S ) 1 (XΦ T S ) 2
problem geometrically amounts to find  corresponding to an unmixing lying within the smallest radius 1 ball. Due to the morphological diversity assumption, it is of uttermost importance to note that such a ball is mainly determined by the highest amplitude source samples. This is the core idea behind GMCA : emphasizing these few high amplitude samples enables an enhanced separation quality. Placing such an emphasis can be done by discarding all the lowest amplitude samples and setting them to zero, which corresponds to setting high hyperparameters M S 8 . We would however suffer from a difficult trade-off :

-If too many samples are kept within the estimation process : taking into account the lowest amplitude samples can worsen the separation results since : i) due to the morphological diversity assumption, the lowest samples are the least discriminative (the support can be very joint between the lowest amplitude samples of the sources) ; ii) the lowest amplitude coefficients are highly likely to be mainly due to the noise N.

-If too few samples are kept : taking only the very highest amplitude coefficients is treacherous since i) the morphological diversity assumption can be an approximation, and small partial correlations 9 can occur even for the highest amplitude samples ; ii) taking into account too few samples can create a lack of statistics, making the problem much more ill-posed.

To bypass this difficult trade-off, GMCA proposes an adaptive hyperparameter M S choice, by making them decrease along the iterations. In the first iterations, only the most prominent samples of A (l-1) † X are kept, and the lower amplitude samples are increasingly added within the estimation process to increase the statistics. Since the M S hyperparameters are thus varying along the iterations l, we shall now denote them as M S (l) in the following.

A final threshold choice based on a fixed point interpretation We now have justified the principle of decreasing thresholds through the morphological diversity. Beyond this general decreasing principle, a natural remaining question is : which value to give to the M S (l) in practice ?

In GMCA, such values are chosen according to a fixed point argument. More specifically, since the algorithm is based on pALS, the thresholding applies to the leastsquare estimate of the sources (cf. Sec. C.4.1). If we assume that after many iterations the algorithm has ultimately stabilized on an estimate Â(∞) close to the true mixing matrix A * , the corresponding source update before thresholding is given 8. Since high M S correspond to a high thresholding S M S in the pALS algorithm. 9. That is, samples with multiple active coefficients.

by10 :

S(∞) = Â(∞) † X A * † X = A * † (A * S * + N) = S * + A † * N (III.32)
Where A † * N is a Gaussian noise (since it was assumed that N was Gaussian).

That is, the sources before thresholding S(∞) are equal to the true ones S * up to an additive Gaussian noise. Thus, the M S (∞) hyper-parameter choice can be understood as choosing the hyper-parameters of a sparse signal Gaussian denoising problem, which has been well studied.

More specifically, let us assume that no reweighted 1 is used and that M S

(∞) = Diag(µ (∞) 1 , µ (∞) 2 , ..., µ (∞) n )1 n×t . For each source S(∞) i , i ∈ [1, n],
the threshold µ i is chosen such that it aims at setting to 0 the small coefficients of S(∞) that should mainly correspond to the noise (A † * N) i . Such a choice is usually performed through a detection procedure using the "κσ " rule, κ ∈ R + . For instance, if κ = 3, the probability that a coefficient of S(∞) i with a larger amplitude than 3σ i , with σ i the standard deviation of the Gaussian noise (A * † N) i , corresponds to noise only is roughly 0.4%. A first practical difficulty is however that the standard deviation σ i of the noise (A † * N) i is unknown (if only because A * itself is not known). It can fortunately be approximated using the Median Absolute Deviation, defined as :

∀u ∈ R t , mad(u) = median i∈[1,t] |u i -median i∈[1,t] (u i )| (III.33)
We further extend this definition to matrices, by taking their row-wise mad : MAD :

R n×t → R n , such that ∀i ∈ [1, n], ∀U ∈ R n×t , MAD(U) i = mad(U i ).
In our case, we have that σ i 1.48×MAD((A † * N) i ). Furthermore, since S * is assumed to be sparse, the MAD operator is quite insensitive to it, and it is thus possible to directly estimate κσ i from S(∞) , as κ MAD( S(∞) )

κ MAD(S * + A † * N) κ MAD(A † * N) = (µ (∞) 1 , µ (∞) 2 , ..., µ (∞) n ) T .
A second, more challenging difficulty is that we have made the assumption that the hyperparameter M S (∞) were chosen according to the point where the algorithm stabilizes, which rises two issues :

-This requires to know the result of the algorithm before launching it ; -This further requires that the result of the algorithm will be good in the sense that Â(∞) A * , which is not trivial since A * is precisely unknown. To bypass such pitfalls, the hyperparameter choice in GMCA is done using the current, at iteration l, estimation of S(l) :

µ (l) 1 , µ (l) 2 , ..., µ (l) n T = κ × MAD( Â(l-1) † X) (III.34)

54

Chapitre III. Optimization frameworks for sparse BSS This is particularly appealing since :

-It alleviate the need to know the result of the algorithm ;

-It draws on the conclusions of the previous subsection, namely that an adaptive decreasing threshold choice enables to take into account the morphological diversity assumption. Indeed, we have that :

µ (l) 1 , µ (l) 2 , ..., µ (l) n T = κ×MAD( Â(l-1) † X) = κ×MAD( Â(l-1) † A * S * + Â(l-1) † N) (III.35) Thus, as during the first iterations Â(l-1) † A * = I d , the term MAD( Â(l-1) † A * S * ) is non zero (it is the MAD of a mixing of sparse signal, which is itself non- sparse), implying high µ (l) 1 , µ (l) 2 , ..., µ (l)
n . When the estimation improves, Â(l-1) † A * becomes closer to the identity and thus the term Â(l-1) † A * S * becomes sparser, and the corresponding µ

(l) 1 , µ (l) 2 , ..., µ (l) n lower.
GMCA, summary GMCA [Bobin et al. 2007] can be written in the following way. The structure is similar to pALS, but the scheme is enhanced with an automatic hyperparameter choice :

GMCA( Â(0) , Ŝ(0) )
Requires only X While not stabilized or maximum number of iterations not reached, iterate over (l) :

1 -Automatic parameter choice :

µ (l) 1 , µ (l) 2 , ..., µ (l) n T = κ × MAD( Â(l-1) † X) (III.36) M S (l) = Diag(µ (l) 1 , µ (l) 2 , ..., µ (l) n )1 n×t (III.37)
2 -S is updated using the current Â(l-1) .

Ŝ(l) = S M S (l) Â(l-1) † X (III.38) 3 -A is updated using the current Ŝ(l) : Â(l) = Π . 2 =1 X Ŝ(l) † (III.39)
Lastly, we would like to highlight that other, more advanced scheme for setting the hyperparameters M S of GMCA have been proposed (see Appendix D for a presentation of some of these extensions). We do not however focus on this aspect, since the common rationale is the same : use a decreasing hyperparameter choice and set the final values according to the fixed point argument.

C.5 AMCA, an extension of GMCA

As evoked in the previous section and already relatively visible in Figure III.2, the morphological diversity hypothesis might be a too strong assumption in practical cases, for which some samples can have several high amplitude coefficients. In this situation GMCA, that grounds its hyperparameter estimation on such high amplitude samples, can have lowered performances. The Adaptive Morphological Component Analysis (AMCA - [Bobin et al. 2015]) has been introduced to tackle such difficult cases. The main idea is the following : AMCA assigns a weight to each sample X j , j ∈ [1, t] which will enable to discard in the estimation process samples suffering from partial correlations : that is, the ones that do not respect the morphological diversity. AMCA weights are grouped as the diagonal elements of a diagonal matrix M ∈ R t×t . Let us write Ω C the source samples for which the morphological diversity does not hold, and Ω * the samples for which it does holds. If such sets are known beforehand, creating the M matrix is easy : for samples X j ∈ Ω C , M jj = 0 (thus discarding them) and for samples X j ∈ Ω * , M jj = 1 (thus keeping them unchanged). Instead of minimizing Eq. II.8, the new cost function 11 is then given by 12 : argmin

A∈R m×n ,S∈R n×t 1 2 Tr (X -AS)Φ T S MΦ S (X -AS) + R S (SΦ T S ) 1 +ι {∀i∈[1,n]; A i 2 2 =1} (A) (III.40)
Where, compared to Eq. II.8, the data fidelity term has been changed to take into account the weights M. Since Ω C and Ω * are however unknown in practice, the authors proposed to relax this method using an adaptive reweighting procedure. They thus use a M matrix chosen based on the sparsity level of SΦ T S columns :

M = Diag j∈[1,t]       1 (SΦ T S ) j p (l) (SΦ T S ) j 2 + ε       (III.41)
where ε is a small constant that may be required when some columns have vanishing p norms, and p ∈ [0, 1] is a parameter determining the influence in the reweighting of the the source samples sparsity level. In brief, p changes during the iterations l : it is close to 1 at the beginning, to avoid to penalize too strongly imperfectly unmixed sources, and it decreases during the following iterations when the unmixing becomes better.

Chapitre IV

Using PALM in Sparse BSS

A Introduction

As mentioned in the previous chapter, the optimization strategy is of uttermost importance with regards to the performances of sparse BSS. Among the algorithms presented earlier, the GMCA algorithm [Bobin et al. 2007] (cf. Sec. III-C.4.2) has well established over the last decade its ability to handle various problems. This in particular due to its heuristic parameter choice -described in Sec. III-C.4.2 -making them decrease during the iterative process and benefitting from the morphological diversity assumption. The goal of this chapter is however to depart from this approach and give a closer look to algorithms that have been introduced more recently and that are usually used with fixed parameters. More specifically, we will focus on PALM ( [START_REF] Bolte | [END_REF] -cf. Sec III-C.3.3), both exploring empirically its performances in the context of sparse BSS, and trying to draw the quintessence from it by re-using heuristics in the same spirit as in GMCA.

A.1 What is a good sparse BSS algorithm ?

But first, a questions is : what is a "good" sparse BSS algorithm ? Generally speaking, any sparse BSS algorithm yields an estimate Θ of the true mixture parameters Θ * = {A * , S * }. This estimate depends on the data X : of course, through Θ * but also through the noise N. To find a good estimate Θ, a practitioner has access to a given data set X but can choose a solver A, an initial point Θ(0) and a set of regularization parameters R S , which can be formulated as follows :

Θ = A Θ(0) , R S ; X
The goal of this chapter is to introduce a reliable, effective and versatile algorithmic framework to tackle sparse BSS problems, which is needed for real-world large-scale applications :

-Reliability means that the solution has a low dependence on the initial point, which is important when little is known about the solution. More specifically with given R S and X, the variance over Θ(0) of the separation quality (measured by an estimator C A -see Appendix B for the one used here) is as small as possible.

-Efficiency has to do with the algorithmic framework and the regularization parameter tuning strategy : for a given dataset X, it should be possible to obtain a solution close to the true Θ without having to perform a cumbersome regularization parameter choosing. In particular, this is achieved if the algorithmic framework is not too sensitive to the regularizing parameter choice (the variance of C A over R S is small for a given dataset X and initialization Θ(0) ) or if some good automatic thresholding strategy is available.

-Versatibility in sparse BSS deals with the fact that the efficiency must hold for various data X settings : the thresholding strategy must generalize well to different data. In the case of a threshold choice performed through a mere grid search, versatility can be assessed by looking at the variance of C A over different datasets X for given Θ(0) and R S . If such a variance is low, it means that a threshold that has been determined for a specific dataset can be re-used easily for another one.

A.2 Questions arising from PALM use in sparse BSS

While the PALM algorithm has become one of the most attractive optimization frameworks for tackling generic matrix factorization [START_REF] Repetti | [END_REF], Lanaras et al. 2015, Thouvenin et al. 2015, Bao et al. 2016, Pierre et al. 2015], its application to sparse BSS raises open questions :

-What are the major limitations of minimizing Eq. (II.8) using the PALM algorithm to perform sparse BSS ? The investigations presented in the following highlight the low reliability of PALM, while performing sparse BSS by minimizing Eq. (II.8) using it exhibits a dramatic lack of efficiency and versatility. The outputs of Eq. (II.8) are indeed highly sensitive to the regularization parameters in some scenarios. In these, PALM -and potentially any other algorithm -is therefore impractical if not associated with an automatic regularization parameter strategy. This is all the more problematic as the regularizing parameters choice is not always well discussed in the literature (although it is a difficult problem [Mensch et al. 2018], in many works it is either not discussed, or sums up to a grid search, or uses the true factorization [START_REF] Repetti | [END_REF], Lanaras et al. 2015, Thouvenin et al. 2015, Bao et al. 2016]).

It is however also empirically shown here that potentially accurate results could be found if optimal regularization parameters for a given (good) initialization were known beforehand.

-Is it possible to improve the estimation yielded by the minimization of Eq. (1) using PALM with heuristic techniques ? The hope is then twofold : i) reach PALM potential accurate results when adapted regularizing parameters are used ; ii) benefit from PALM mathematical guarantees. It has to be emphasized that using heuristics in other matrix factorization algorithms has already known a wide success, both for sparse BSS [Bobin et al. 2007] and for NMF [Vandaele et al. 2016]. We show how combining PALM with heuristic techniques can enhance the separation quality by providing adapted initializations and regularization parameters.

B Chapter outline

In section C, we investigate the practical limitations of performing sparse BSS by minimizing Eq. (II.8) using PALM. In section D, we show and explain the difficulty of straightforwardly extending existing heuristic techniques to PALM to make it easier to use. In section E, we re-use the results of section C and D to handle each element hindering PALM applicability in sparse BSS and, building on them, to rationalize a 2-step approach (with an initialization stage followed by a refinement procedure). In section F, the quality of the 2-step approach is demonstrated on realistic astrophysical data, while in section G, the limitations of the approach are studied. Several extensions of this work are furthermore shortly discussed in Appendix E (in particular, a quick study of BCD is performed, and other 2-step approaches are proposed).

C Limitations of minimizing Eq. (II.8) with PALM to perform sparse BSS, an empirical study

The objective of this section is to empirically shed light on the limitations of using the PALM algorithm along with Eq. (II.8) to handle the sparse BSS problem. For that purpose, we will evaluate the accuracy of the final point estimate Θ = { Â, Ŝ} yielded by PALM in terms of an estimator (which can be computed here only due to the fact that the true matrices are simulated and therefore known). More specifically, we will highlight the limitations of using PALM to minimize problem (II.8) to perform sparse BSS in terms of efficiency, reliability and versatility.

At this point, it might be important to highlight two elements :

-We study the results in terms of separation quality measured by an estimator C A , and not only in terms of the cost function (II.8). That is, we try to find a good critical point corresponding to a true physical factorization.

-The phenomena we study here could be inherent to the use of Eq. (II.8) for performing BSS, since the role of PALM is only to perform its minimization. However, in the light of the previous remark and since (II.8) is non-convex, the use of another minimization scheme could lead to a different (local) minimum with a different quality in terms of the estimator. For instance, another minimization algorithm could be more or less sensitive to the regularization parameter values 1 . In this work, we mainly study the separation quality when minimizing Eq. (II.8) with PALM (and therefore not problem (II.8) in general -however, preliminary results with BCD are also presented in Appendix E, tending to suggest that several conclusions might be quite inherent to the cost function).

C.1 Gist of the experiments : what is to be studied ?

The goal of this subsection is to determine which factors have a potential influence on the final separation quality and must therefore be taken into account in the study of PALM for sparse BSS. This can be inferred by detailing the first update step of the sources in the PALM algorithm :

S = Ŝ(0) - γ Â(0)T Â(0) 2 Â(0)T ( Â(0) Ŝ(0) -X) = Ŝ(0) + γ Â(0)T Â(0) 2 ( Â(0)T A * S * -Â(0)T Â(0) Ŝ(0) + Â(0)T N) (IV.1)
The gradient descent step is followed by the application of the proximal operator of the penalization term :

S (1) = S γR S Â(0)T Â(0) 2 S (IV.2)
Therefore, the best estimation of S * will potentially depend on the starting point Θ(0) = { Â(0) , Ŝ(0) }, the true mixture parameters Θ = {A * , S * }, the noise N and the regularizing parameters R S .

More specifically, to relate this discussion to the desired properties of a sparse BSS algorithm and outline the experimental protocol :

-The efficiency will be studied by trying several values of R S for a given experimental setting. This will give an insight of the estimate quality sensibility (in terms of C A ) yielded by the minimization of Eq. (II.8) using PALM.

-The reliability will be studied varying the initial points Θ(0) and more specifically Â(0) .

-The versatility will be studied by changing the experimental setting : Θ = {A * , S * }, and the noise N. More specifically, 4 types/cases of experiments will be described in Sec. C.2 and we will analyse :

-The inter-case versatility : for instance, how a change in the source distribution affects the estimate quality ?

-The intra-case versatility : even focusing on a specific case, does using another random realization of Θ = {A * , S * }, and N can affect the quality ?

C. Limitations of minimizing Eq. (II.8) with PALM to perform sparse BSS, an empirical study 61

C.2 Description of the data

To bring out the mechanisms at stake, we will perform experiments using 4 different datasets. In each case, there is n = 2 sources2 .

-Case 1 : The sources are assumed to be exactly sparse in the direct domain and to follow a Bernoulli-Gaussian distribution : a proportion p = 0.1 of the t = 500 samples is non-zero and drawn according to a standard normal distribution. The sources are equilibrated, i.e. they have equal variances. Their dynamic (maximum minus minimum value) is circa 0.6. This example therefore corresponds to very simple synthetic data. The mixing is performed through a matrix A * drawn randomly following a standard normal distribution and modified to have unit columns. Its condition number is C d = 10 and we focus on the exactly determined case : there is an equal number of observations and sources m = n. To complete the creation of the X data, a Gaussian noise is added to the mixing, such that the Signal to Noise Ratio is SNR = 60 dB.

-Case 2 : The sources are assumed to be approximately sparse in the direct domain and to follow a generalized Gaussian distribution of parameter α = 0.25. There is again t = 500 samples and the sources are still equilibrated with a dynamic of circa 0.6. The mixing is performed with the same A * matrix as in case 1 and an additive noise N is added in a similar way. Compared to case 1, case 2 corresponds to a more realistic setting, since the wavelet coefficients of natural images would have a similar distribution.

-Case 3 : The sources are constructed in a similar way as in case 2. However, the noise energy of the first observation is twice the one of the second observation. Furthermore, A * is taken orthogonal (that is, A * T A * = A * A * T = Id) to ensure that the noise projection on the source space has different amplitudes. While this setting might be simpler than the one of cases 2, it allows to study the impact of the noise on the parameter choice.

-Case 4 : The sources come from simulations obtained from real data of Cassiopeia A supernova remnant. These data originate from the Chandra3 X-ray observatory. The sources in these wavelength values correspond to the thermal emission and the iron 2. As displayed in Fig. II.3, each of them consists in a 2D image of resolution t = 128 × 128 pixels, supposed to be approximately sparse in the starlet domain (in which their dynamic range is circa 4 × 10 -3 and 6 × 10 -4 respectively). The mixing is performed with the same A * matrix as in case 1. A Gaussian additive noise is added to the mixing, such that the SNR = 30 dB. Beyond being still more realistic than the other cases, case 4 involves non-equilibrated sources.

C.3 Evaluation protocol

For each experimental scenario, a large number of values for the regularization parameters R S are tested. More specifically, since no reweighting 1 in problem (II.8) is used in this part, the work will be done on Λ S = R S . As we will focus on the case of n = 2 sources, there will be 2 parameters to test, λ 1 and λ 2 . For each (λ 1 , λ 2 ) value, a criterion measuring the separation quality is evaluated. Since a change in the regularization parameters directly impacts the source estimation, we will use the mixing matrix criterion described in [Bobin et al. 2015] (cf. Appendix B) :

C A = mean(|P Â † A * -I d |) (IV.3)
With A * the true mixing matrix and P Â † the pseudo-inverse of the estimate corrected through P for the scale and permutation indeterminacies. The mean is the average of all the elements inside the matrix. For the sake of clarity, the plots will display -10 log 10 (C A ). As such, the higher the values, the better the separation.

C.4 Results and interpretation

C.4.1 Efficiency : sensitivity to the regularization parameters

In this subsection, we study the efficiency, that is we want to get an idea of the sensitivity for a given X of PALM estimate to the regularizing parameters Λ S = R S . Therefore, we take one fixed experiment of each case and we try several (λ 1 , λ 2 ) values. For each of them, the mixing matrix criterion C A is computed from the estimate and reported in a 2D plot, displayed in Figure IV.1. To try to separate as much as possible the efficiency issue to the reliability issue, we launched the algorithm with 5 different random initializations Â(0) and the median is displayed. The first observation corresponds to the high variations induced by changes of regularizing parameters, with a dynamic range of more than 40 dB for case 1, 30 dB for case 2 and 20 dB for case 4. The results change from good to catastrophic, making the choice of the regularizing parameters a crucial point, especially as these high variations are very fast on the plot. In case 1, since the sources have identical distributions and the noise is supposed to be white, it would be expected from the best regularization parameters λ 1 and λ 2 to have similar values. However, while staying close to the diagonal could restrict the space of possible good parameters, it is unfortunately dangerous to fully restrict the search to the diagonal, as demonstrated in Fig. IV.2. In case 1, just using the parameters on the diagonal closest to the best ones of Fig. IV.2a yields a drop of more than 10 dB (which might be linked to a too small number of samples for the noise to be perfectly white). In the non-equilibrated case 4, while not very visible in Even if we assume that it is possible to restrict the regularizing parameter search to the diagonal, the variations are still quick in comparison to the dynamic of the sources. In the simple case 1, a shift of 0.1% of the dynamic of the sources in the regularizing parameters in a parallel direction to the diagonal yields a 7 dB loss. In case 2, a shift of 5% yields a 8 dB drop and for case 3, a 1% change implies a 8 dB drop. This highlights that, for a single experiment (i.e. for a given A * , S * and N), the separation performances are highly sensitive to the choice of the regularization parameters.

C.4.2 Versatility : impact on the regularization parameters of A * , S * and N

Inter-case versatility

We first look at the main potential source of lack of versatility, that is when the distributions of the matrices A * , S * and N change. To do that and for computational reasons, we will merely compare the results displayed in -Impact of S *

The shapes of the plots of the 4 different cases strongly differ. In particular, the observations we made in the previous subsection show that the best regularizing parameter values are highly dependent on the true source S * distributions :

C. Limitations of minimizing Eq. (II.8) with PALM to perform sparse BSS, an empirical study 65

-In addition to remove the noise, the regularizing parameters should limit the presence of remixing or interferences between the sources, which depend on the distribution of S * , the mixing matrix A * as well as the initialization.

-The resulting thresholding induces a bias, called artifacts, which also depend on the distribution of S * . These are also problematic as they can eventually convert into interferences between the sources in the iterative optimization process.

-Impact of the back-projection of N At the vicinity of the true mixture parameters A * and S * , the only source of error that contributes to the estimate of S * originates from propagated noise. The role of the regularization parameters is then merely to avoid a deterioration of the estimate by the noise, that is to threshold the update

A * T (X -A * S * ) = A * T N.
Consequently, an intuitive choice for the regularization parameters for each sources i would be λ i = (A * T N) i ∞ , which implies a clear dependency between the thresholds and the back-projection of N through A * T . This can be studied through case 3 : in this setting, A * is orthogonal and the noise N energy is different over the observations. The thresholds enabling the convergence should therefore be different for the two sources : λ 1 = λ 2 . This is confirmed by Figure IV.1c, which shows that contrary to the other cases, the highest separation quality lies far from the diagonal. Thus, unbalanced backprojected observation noise (that is, with different variances on each source) can make the threshold choice much more difficult4 than in settings where the noise is equilibrated due to optimal parameters further from the diagonal.

Intra-case versatility

We now look at the intra-case versatility, that is we study the impact of changing the realization of the data for a fixed case -more specifically, we will focus on case 1 and 2. To do that, we draw for both cases 10 new random A * , S * and N matrices and create new data. For each of these data and to try to separate the versatility issue from the reliability one, we try 10 initializations and take the median over them. Therefore, we get 10 plots similar to the ones of In the exactly sparse case 1, the dynamic of circa 40 dB is huge compared to the results yielded by the best value of some settings (which are in some realizations of 24 dB). For the approximately sparse case 2, the dynamic is relatively similar. It makes that even without changing the distribution of the random matrices (said differently, even if we know their distributions), choosing good (λ 1 , λ 2 ) values is very difficult because this choice is extremely dependent on the specific data realization that must be handled.

Consequently, the best regularizing parameters are dramatically dependent on the data X, which further highlights the low versatility of performing sparse BSS by minimizing Eq. (II.8) using PALM without any automatic regularization parameter choice. In particular, the extra-case versatility has pointed out that the regularization parameters should be chosen based on S * as well as the backprojected noise.

C.4.3 Reliability : impact of the initialization Â(0) , Ŝ(0)

The impact of the initialization has two theoretical groundings :

-As problem (II.8) is not convex but multi-convex, the algorithm performing its minimization can be trapped in spurious critical points depending on the initial matrices Â(0) and Ŝ(0) .

-The initialization directly impacts the quality of the estimate yielded by specific thresholds through Eq. (IV.2) and the threshold choice. Said differently, for a given initialization, it might be desirable to change the cost function via the choice of thresholds to avoid a spurious critical point that would not be an issue for another initialization.

To quantify this impact, the dynamic of C A over 5 different initializations Â(0) is plotted for case 4 in This is mainly due to regularizing parameters outside the diagonal, which unfortunately correspond to a high separation quality for some experiments. This means than choosing such parameters in these experiments can lead to almost the best performances, as well as very bad ones (3 dB), making the regularization parameter choice still more difficult. These experiments emphasize that the quality of PALM results is highly connected to the initialization. More strikingly, the choice of good regularization parameters in terms of C A is also connected to the initial point. In practice, to avoid relaunching the algorithm with different initializations for given regularization parameter values (which is especially important with large-scale data), it is important to make PALM estimate, as well as the regularization parameter choice, more robust to the initial conditions.

D Enhancing PALM with heuristic techniques

To bypass the cumbersome parameter choice described in the previous section, the goal of this section is to try and show the issues of a straight adaptation within PALM of the heuristic parameter choice used in GMCA. In the first subsection, GMCA heuristic is reminded, and some general remarks about the performances of GMCA made. In section D.2, the heuristic parameter choice is adapted to the PALM framework. In section D.3, the heuristic is tested and shown to work badly within PALM, while in D.4 these empirical results are explained. Tableau IV.1 -GMCA results in terms of C A .

D.1 GMCA heuristic and performances

D.1.1 Reminder concerning GMCA heuristic

From Chapter III, we would like to recall that GMCA proposes an automatic threshold choice by using the following κ-MAD rule :

M S = Diag µ (l) 1 , µ (l) 2 , ..., µ (l) n = Diag(κ × MAD( Â(l-1) † X)) (IV.4) D.1.2 Remark about the performances of GMCA
The practical success of the GMCA algorithm relies on its good reliability as well as the heuristic to automatically tune M S , enabling both efficiency and versatility.

To give an idea of the performances of GMCA in the experimental setting of Sec. C, we launched it on cases 1 and 2, for 10 different initializations and 10 A * , S * and N settings. The results are summarized in Table IV.1.

Compared to PALM, the dynamic over different realizations of A * , S * and N is much lower than the one of PALM, and the worst values are still good, meaning a much higher versatility. The automatic parameter choice yields a high mean of C A , showing GMCA efficiency. The reliability seems to be extremely high as shown by the very low standard deviation over the initializations. However, the results of GMCA are not always as good as the ones yielded by PALM with the best parameters. This might be due that the use of the pseudo inverse A (k) † in GMCA can strongly increase the noise back-projection when Â(l) is badly conditioned. As an example of such a difference between GMCA and PALM, the results of GMCA are 36.3 dB in the experience of case 1 displayed in Fig. IV.2, to be compared with a little more than 40 dB for PALM. This partially justifies the will to make PALM easy to use : to benefit from its potential high accuracy when it is well tuned.

D.2 Adaptation of GMCA heuristic parameter choice to PALM : a first idea

The goal is here to explain a straightforward adaptation to PALM of GMCA automatic parameter choice. For the sake of simplicity, let us assume that there is no reweighting : R S = Λ S . The update of S by PALM in Eq. (III.25) is then :

Ŝ(l) = S γΛ S Â(l-1)T Â(l-1) 2 SPALM (IV.5)
Under the same assumptions, the update of GMCA is the following (with SGMCA = Â(l-1) † X) :

Ŝ(l) = S M S SGMCA (IV.6)
Taking into account the differences between SPALM and SGMCA , the parallel between Eq. (IV.5) and (IV.6) suggests that Λ S could be choosen similarly as the parameters in GMCA (in which case the parameters change over the iterations, but are fixed at the end to ensure the convergence of the algorithm) :

γ A * T A * 2 (λ 1 , λ 2 , ..., λ n ) T = κ × MAD SPALM (IV.7)
If so, and if PALM has converged to both the true matrices A * and S * , then :

γ A * T A * 2 (λ 1 , λ 2 , ..., λ n ) T = κ × MAD S * - γ A * T A * 2 A * T (A * S * -X) κ × γ A * T A * 2 MAD(A * T N) (IV.8)
Where the last line is obtained because S * is assumed sparse and the MAD is robust to outliers. Therefore, using the MAD enables a thresholding of a projection of the noise N, which yields a similar interpretation as in GMCA (cf. Sec. III-C.4.2). This parallel must however be tempered :

-It only holds when and if PALM has converged towards good A and S ; -The projection is not performed through the pseudo-inverse as in GMCA, but rather through A * T . Both projections however merge when A is orthogonal.

Remark : Choosing the regularization parameters in general inverse problems has been the subject of several studies, which propose various ways for setting them. Among them, one can highlight the Stein's Unbiased Risk Estimator -SURE -method (and its extended versions [Stein 1981, Eldar 2009, Giryes et al. 2011, Deledalle et al. 2014]), the generalized cross-validation [Golub et al. 1979[START_REF] Lukas | [END_REF]], the L-curve [START_REF] Hansen | The use of the L-curve in the regularization of discrete ill-posed problems[END_REF], the discrepancy principle (or some variants [Almeida & Figueiredo 2013]) or some Bayesian methods [Pereyra et al. 2015, Vidal & Pereyra 2018]. However, most of these are not directly tractable in the case of sparse BSS for several reasons : i) we are dealing with the blind setting, in which the linear operator A * must also be evaluated, making the use of these methods more complicated ; ii) some of these methods can be computationally expensive in the large-scale setting, since they require to try several regularization parameters and compute a criterion to decide which one to choose ; iii) some of them, such as SURE, use as a criterion for the parameter choice the estimated MSE : it is not clear if this criterion is the most relevant one in sparse BSS (cf. [START_REF] Feng | Revisiting sparse ICA from a synthesis point of view : Blind source separation for over and underdetermined mixtures[END_REF]) ; iv) some of them have mostly been applied in the case of linear solutions (e.g. by using a 2 regularization instead of 1 -this has been emphasized in the context of deconvolution in [Almeida & Figueiredo 2013]). Therefore, choosing the regularization parameters using the MAD seems one of the most interesting solutions since it has already lead to enhanced separation quality within the context of GMCA. In Appendix E, we however try one of the previous methods in sparse BSS to back this claim.

D.3 Illustration

The automatic parameter setting based on the MAD described in the previous subsection is tried inside of a PALM algorithm for different κ values. The parameter κ is the same for both sources, since the MAD is supposed to be adaptive enough to the noise. The metric on the separation quality C A is the same as used in Sec. C and the data is the same as in case 1 : it is obtained from two exactly sparse sources with a square mixing matrix having a condition number of 10.

D.4 Understanding the limitations of the heuristic in the scope of PALM

We previously pointed out that the thresholds are not identical in GMCA and PALM. A key difference is that the thresholding applies to different quantities :

-In GMCA, it is applied on a least-square estimate of the sources after a direct inversion of the mixing matrix ;

-In PALM, it applies on a single gradient descent step, that is during the inversion.

Therefore, the estimation errors originating from both an imperfect unmixing and the noise are different in the two algorithms. Since it is the role of the thresholds to filter out these estimation errors, they do not have the same impact and the optimal ones may differ from an algorithm to the other.

To better understand this role, let us assume that both GMCA and PALM algorithms are initialized with the true mixing matrix A * only. Let us further assume on the contrary that the initialization of the sources is not perfect and can be written as S = S * + s with s the error made on the sources. For GMCA, the thresholds given by the MAD heuristic then are :

(µ 1 , µ 2 , ..., µ n ) T = κ MAD(A * † X) = κ MAD(S * + A * † N) κ MAD(A * † N) (IV.9)
Where the last equality is obtained because the MAD operator is not sensible to sparse signals. In this case, the thresholds are thus set according to the backprojection of the noise only and an imperfect estimation of S * does not affect the thresholding strategy as long as the mixing matrix is well estimated. In contrast, when using PALM the thresholds are given by :

γ A * T A * 2 (λ 1 , λ 2 , ..., λ n ) T = κ × MAD S - γ A * T A * 2 A * T (A * S -X) = κ × MAD S - γ A * T A * 2 A * T (A * (S * + s) -X) κ γ A * T A * 2 MAD A * T N -A * T A * s (IV.10)
The previous equation highlights the following (and unfortunately complementary) issues of using the MAD inside of PALM :

-An inappropriate threshold choice :

Compared to GMCA, the thresholds are calculated with the additional detrimental A * T A * s interference term. Indeed, they are computed after only one gradient step and therefore when the interferences are still high. Another interpretation is that the MAD is computed on an approximation of the minimization of the data fidelity term, namely a gradient step update, and its results are therefore not as accurate as in GMCA.

This hinders the interpretation of the use of the MAD in terms of noise removal, particularly when the remixing A * T A * s has more energy than the noise term A * T N. More specifically, this is the case when the columns of A * are strongly correlated, as testified by the product A * T A * , or when the error on the sources s are high (which is in particular the case if the initialization is random or if there are many artifacts). If s is further strongly non-Gaussian, the MAD use is still less relevant (at least, when the number of sources n is low).

-Higher interferences in the estimation process :

The gradient update on the sources creates at each iteration a remixing of the error s due to the A * T A * s term. While this could be a limited issue if A was fixed, since in PALM only one proximal gradient step is performed on S before the update of the mixing matrix, A will be computed from sources with high interferences.

In particular, this issue is very important with high thresholds that induce artefacts and therefore a large s, creating more interferences (which in turn feeds the problem of the previous point, creating more artefacts and so on...). It makes that high thresholds are not that valuable in PALM. Therefore, contrary to GMCA in which the A * T A * s is not present, it is harder to exploit the Morphological Diversity assumption. Indeed, while the most significant coefficients are still assumed to be the most interesting for the separation, selecting only them will create interferences counterbalancing their positive effect.

To conclude, in PALM high thresholds are less interesting than in GMCA and using the MAD heuristic does not results in good separations in practice if the level of the interferences is not negligible in comparison to the noise contribution.

E Combining GMCA and PALM : a hybrid strategy

The previous part has demonstrated that the introduction of an heuristic based only on the back-projection of the noise does not yield satisfactory results within a single PALM. We therefore rationalize an approach tackling the various elements (cf. Sec. C) hindering the applicability of PALM in sparse BSS by combining the best of GMCA and PALM in a two-step approach. The first warm-up stage consists in a GMCA. It is followed by a refinement stage during which PALM is performed retaining as much information as possible coming from the warm-up stage. Inside the refinement stage, the MAD heuristic is used to choose PALM parameters. The main idea of the two step approach is to reduce the interferences by a good initialization of the refinement stage and the artifacts using a reweighting information coming from the warm-up stage. It has to be emphasized that this 2-step approach has already been empirically used in two previous works [Chenot & Bobin 2018,Kervazo et al. 2018]. Therefore, the novelty of the following is to provide an in-depth justification of this approach while the previous results were only empirical.

More specifically, this justification describes in details all the required elements for this approach to work, which is done both with theoretical arguments along with experiments to confirm these on new datasets.

E.1 Motivation of the two step approach

Our 2-step approach has several motivations :

-Dealing with the factors deteriorating the results of sparse BSS (in terms of versatility, efficiency and reliability) when using PALM for minimizing Eq. (II.8) :

While the MAD heuristic enables the refinement stage to handle the noise back-projection, using some reweighting information from the warm-up stage enables to take into account the distribution of S * , which was the second factor of the lack of versatility of using PALM in the context of sparse BSS (cf. Sec. C). The automatic parameter choice then enables to circumvent in the refinement stage the difficult parameter choice and the lack of efficiency. Finally, using a reliable algorithm such as GMCA as a warm-up stage enables a reliable global 2-step algorithm. Furthermore, it enables a good initialization of PALM (which is reached using decreasing regularization parameters and benefiting from the morphological diversity, in contrast to what is usually done when using a single PALM with fixed regularization parameters).

-Enabling to use the MAD heuristic : One of the main advantages of combining GMCA and PALM is to initialize PALM with the output ÂGMCA and ŜGMCA of GMCA. Since both ÂGMCA and ŜGMCA are close to A * and S * , the level of the interferences is relatively low in comparison to the noise. Therefore, following the conclusion of Sec. D.4, the MAD is still relatively accurate to derive the thresholds once at the beginning of the refinement step and the interpretation in terms of noise removal is more accurate than with a poor initialization, due to smaller interferences. In addition, the interferences are further indirectly reduced by the use of the reweighting, which reduces the artefacts that would be transformed into interferences by the gradient step.

-Keeping mathematical guarantees and PALM high potential accuracy : While GMCA is only a proxy, the 2-step algorithm will attempt to solve exactly Eq. (II.8) as PALM does. It will further benefit from PALM potential high accuracy (cf. Sec. D.1.2). Moreover, since PALM is proved to converge under mild assumptions [START_REF] Bolte | [END_REF], so is the 2-step algorithm 5 .

5. Note that here only a single PALM is used in the refinement step. It could be interesting to study a two-loop algorithm, in which several PALM would be launched until stabilization (if any), with each time new regularization parameters computed on the solution of the previous PALM. It seems however clear that such a process would preclude any convergence guarantee of the whole algorithm, as it is not straightforward that the regularization parameters would converge.

E.2 Use information on S * from GMCA : reweighted 1 As evoked in the previous subsection, it is possible to take into account in the refinement stage the information from S * coming from the warm-up stage. To that end, minimizing the artifacts can be carried out by improving the sparse regularizer. This is done by resorting to a reweighted 1 regularizer [START_REF] Candes | [END_REF].

In this setting, one can benefit from the first guess estimate ŜGMCA to compute the reweighting matrix G in problem (II.8). This generally leads to a significant decrease of the artifacts, which eventually reduces the importance of the interferences introduced by the A * T A * s term. This will improve the efficiency of the proposed hybrid heuristic, which interpretation is based on noise removal only. In this work, we use the following reweighting scheme :

G j i = ε ε + | Ŝj GMCA i | ŜGMCA i ∞ (IV.11)
where ε is a small constant (here 10 -3 ), G j i the coefficient of G corresponding to the i th line and j th column and Ŝi the i th line of Ŝ. In brief, the thresholds are lowered for the largest samples of the estimated sources, reducing the bias introduced by the soft-thresholding and therefore the error s and the interferences.

E.3 Final algorithm

The previous remarks lead to the following algorithm using a smart initialization and thresholding strategy : Input : X (data matrix) -Warm-up stage :

Random initialization Â(0) and Ŝ(0) ÂGMCA , ŜGMCA = GMCA(X, Â(0) , Ŝ(0) )

-Refinement stage : Update of the reweighting information in Eq. (II.8) :

G j i = ε ε + | Ŝj GMCA i | ŜGMCA i ∞ (IV.12)
Update of the parameters Λ S in Eq. (II.8) :

γ ÂT Â 2 (λ 1 , λ 2 , ..., λ n ) T = κ × MAD  Ŝ GMCA - γ ÂGMCA T ÂGMCA 2 ÂT GMCA ( ÂGMCA ŜGMCA -X)   (IV.13) PALM step, initialization coming from GMCA : ÂPALM , ŜPALM = PALM(X, ÂGMCA , ŜGMCA )

E.4 Complexity of the algorithm

We here derive the complexity of one iteration of each of the 2 steps of the algorithm.

E.4.1 Initialization stage

Each iteration of the warm-up stage can be decomposed into the following elementary steps : i) the pseudo-inverse is performed using the singular value decomposition of a n × n matrix, which yield an overall complexity of O(n 3 + m 2 n + nmt) ; ii) the thresholding-strategy first requires the evaluation of the threshold values, which has a complexity of nt ; iii) the soft-thresholding step itself has a complexity of O(nt) ; and iv) updating A is finally performed using a conjugate gradient algorithm, whose complexity is known to depend on the number of non-zero entries in S and on the condition of this matrix C d (S). An upperbound for this complexity is O(nt C d (S)). The final estimate of the complexity of a single iteration is thus given by :

m 2 n + n 3 + nmt + nt C d (S) (IV.14)
With C d (S) the condition number of S.

E.4.2 Refinement stage

For the refinement stage, i) the update of S is dominated by the multiplication required for ÂT (X -ÂŜ ), which has a O(nmt) complexity, and the computation of ÂT Â 2 having an overall complexity of O(n 2 m) (for instance, using the power method). The thresholding and weight computation complexities are negligible ; ii) similarly, for A update the complexity is O(n 2 t + nmt). Therefore, the final estimate of the complexity of a single iteration of the refinement stage is given by : nmt

+ n 2 (m + t) (IV.15)
Beyond this complexity, the overall number of required iterations needed to give "good" results for GMCA and convergence for PALM is also of uttermost importance. In practice, the number of PALM iterations tends to be much higher than the one of GMCA. Therefore, the refinement stage is in practice much more computationally expensive than the warm-up one. The whole 2-step algorithm is however much faster than a single isolated PALM. This is due to the warm-up stage that enables the refinement stage to start from a good initialization and therefore to converge in less iterations.

E.5 Illustration

The experimental protocol is the same as described in Sec. D.3 (two exactly sparse sources with a square mixing matrix having a condition number of 10 and a SNR of 60 dB) except that the 2-step algorithm is used instead of PALM. The results are plotted in Fig. IV.5. With values of C A higher than 33 dB, the demixing is close to the best ones obtained with the exhaustive search in Sec. C. Compared to PALM only, the variance of results over different initializations is also largely improved as it is close to zero, which shows the increased reliability of the algorithm with regards to the initialization.

F Application to a realistic data separation problem in astrophysics

F.1 Description of the data

The following numerical experiments are carried out on simulated astrophysical data, which have been generated from real Chandra 6 observations of the Cassiopeia A supernova remnants (cf. Sec. II-A.1.3). The observations are made of a linear combination of three components : the synchrotron emission and 2 redshifted iron (Fe) emission lines 7 . Compared to case 4, there are now n = 3 sources (which creates 6. http ://chandra.harvard.edu/ 7. More precisely, the emissions lines correspond to the telescope impulse response to Diracshaped emission lines, which are modeled as Gaussian-shaped spectra. These lines are centered about different energy values, which depend on the relative speed of propagation of each iron

F. Application to a realistic data separation problem in astrophysics 77

Wavelength value (arbitrary units) Energy (arbitrary units)

F.2 Results

Results for different SNR values in terms of C A are displayed in Table IV.2. The original sources and the ones estimated for a SNR of 30 dB, as well as the difference between both rows, are shown in Fig. IV.7. The "PALM" line corresponds to a PALM algorithm equipped with the MAD heuristic described in Section D (which is already a potential improvement compared to an exhaustive search on the regularizing parameters -in the sense that it makes possible to use PALM in the large-scale context -if the interferences are low compared to the noise level). To make the comparison fair with RNA [Zibulevsky 2003] and EFICA [Koldovsky et al. 2006], in this subsection the data were pre-whitened. The separation quality of the 2-step algorithm is good, both for the estimation of A * (cf. Table IV.2) and S * (cf. Fig IV .7). In this experiment, the 2-step algorithm always obtains better results than both PALM and GMCA with a gain of about 2 dB for all tested SNR. This highlights the advantage provided by the 2-step approach compared to either GMCA or PALM alone. Interestingly, the PALM algorithm alone provides rather reasonable separation recomponent due to the Doppler effect. The synchrotron component has a power emission law. These components are representative of typical supernovae remnants in the energy band 5000 -6000 eV (electron-volt). The standard sparse ICA-based methods, namely RNA and EFICA, perform rather poorly at low SNR, which highlights their higher sensitivity to noise. In contrast, they perform similarly (or even better) to GMCA in the low noise regime. G Discussion on the 2-step approach G.1 Limitations of the current 2-step strategy

While the 2-step approach has been shown to perform very well in a realistic astrophysical setting, its limitations and applicability are now discussed in more details.

G.1.1 Explanation of the limitations on the realistic setting

To further understand in which regimes the 2-step is valuable or not, let us have a closer look at Eq. (IV.10) where two important terms must be highlighted : i) the back-projected error on the first guess sources A * T A * s and ii) the noise backprojection A * T N. For a fixed starting point (i.e. fixed error s), both the condition number of the mixing matrix A * T and the SNR define different regimes. As an illustration, experiments are performed with the astrophysical data described before but with random non-negative mixing matrices with a larger condition number equal to 10. We carry out 10 Monte-Carlo simulations over A * , with for each 10 different random initializations. The results are displayed in Fig. IV.8a, in which the average over the initializations is used. The error bars correspond to the quartiles of the results for different A * . For larger condition numbers, the re-mixing effect in the PALM iterations plays a prominent role since it tends to concentrate the first guess errors in the subspace spanned by the eigenvectors that are associated with the largest eigenvalues of the Gram matrix of A * . Consequently, the MAD heuristic is more likely to yield badly estimated regularization parameters and the results of the 2-step approach compared to GMCA only are worse than in the previous subsection. This gain will likely depend on the relative levels of A * T A * s and A * T N, which are displayed in Fig. IV.8b. When the noise level is large enough, the term A * T N dominates and the proposed MAD heuristic will perform correctly, which will be favorable for the 2-step approach. This can be observed when the SNR is below 25 dB but not too small : for still smaller SNR, GMCA does not work well, which makes that the 2-step approach has deteriorated performances due to a bad initialization (this is probably due to the pseudo-inverse in GMCA of a badly conditioned A in the presence of strong noise). For larger SNR (in the range 25 -50 dB in A * T A * s becomes dominant. In that case, the 2-step approach might degrade the GMCA solution. Finally for larger SNR, the contribution of the noise to the thresholds is almost zero. While the MAD therefore yields a bad estimation of optimal parameters, the thresholds are however low because of low interferences due to the relatively good initial point proposed by GMCA. It makes that the solution yielded by the refinement step does not deviate significantly from a simple minimization of the data fidelity term. Since in GMCA the data fidelity term is minimized at each iteration, the solution of the warm-up stage is close to one of its minimizers and therefore the solution of the refinement step does not change much from the one of GMCA. Said differently, in this regime even if the thresholds are bad, they are low enough for the refinement stage estimate to stay close from the one yielded by the warm-up stage, which is already quite good because GMCA performs well with low noise8 .

G.1.2 Further understanding of the limitations of the two step strategy

To further understand the limitations of the two step strategy and the mechanisms at stake, we propose to come back to simulated data. Doing so should in particular enable to understand more deeply the influence of the two terms A * T A * s and A * T N by playing on the sparsity level, the noise level and the condition number of A * . More specifically, the experiences are similar to the one of case 2 : the sources follow a generalized Gaussian distribution of parameter α, which will vary in the range between α = 0.05 (very sparse sources) and α = 1 to study the impact of s in A * T A * s. The noise is Gaussian, with a SNR that will vary from 10 dB to 120 dB, enabling to study N in the term A * T N. Finally, we will look at the influence of A * in both terms by trying six condition number C d values, namely 1, 2, 5, 10, 20 and 100. The matrix A * is further drawn randomly from a standard normal distribution and modified to have unit columns. For each parameter value, we launch the 2-step algorithm and compute the mixing matrix criterion C A . Each value of C d yields a 2-D image, displayed in Fig. IV.9. The 6 upper plots correspond the results of the GMCA warm-up stage, and the 6 lower plots to the improvement of the PALM refinement stage over the warm-up stage. The differences between lower and upper parts correspond to two phenomena : i) the different behavior between PALM and the pALS scheme ; ii) the relevance of MAD use within PALM. Generally speaking, the differences between GMCA and the two step strategy are increased by higher condition numbers C d , which was expected since in the case of orthogonal matrices (C d = 1), the pseudo-inverse is equal to the transpose :

A † = A T .
Moreover, the plots of Fig. IV.9 confirm what was observed in the previous subsection with realistic sources : the 2-step strategy improves GMCA results when the noise level is high (left part of the plots), that is when GMCA does not work well. This is due to the fact that i) the pALS scheme of GMCA and in particular the use of the pseudo-inverse amplifies the noise, which is not the case in PALM ; ii) the MAD is relevant, since the Gaussian noise A * T N dominates over the interferences A * T A * s, making the gradient step similar to Gaussian if the estimates are close to the true values. Point ii) can further be made more precise by looking at the shape of the first gradient step within the refinement step (that is, the shape of A T (X -AS) at the first iteration of PALM). Such a shape is displayed in .9. On the diagonal, while the residual are not Gaussian they are not too different either, explaining acceptable results of the 2-step approach (although the residual seems to be further from a Gaussian in the lower left plots of Fig. IV.9, the better results are likely to come from the 9. The term "good" is to be understood in comparison to the output of the GMCA stage, and not necessarily in terms of the absolute performances of the whole 2-step algorithm -which can be quite bad in this difficult setting where the sources are not very sparse and the noise level high.

fact that the MAD is fairly insensitive to the sparse contamination coming from the interferences A * T A * s -such a contamination is however less sparse in the upper right corner). Furthermore, we have already emphasized that the use of reweighted 1 enables to lower the artifacts, that might be transformed into interferences by the gradient step in PALM (cf. Sec. E.2). Looking in Fig. IV.11 at the (ideal) residual at convergence A * T (N -A * s) for highly sparse S * and high SNR enables understanding another interest of reweighting : it transforms a highly non-Gaussian residual into a Gaussian one. Using reweighted 1 , the MAD is therefore more likely to perform correctly and to enables a threshold choice cancelling the erroneous updates, which confirms that the reweighting is of uttermost importance in the 2-step strategy.

Conclusion

The applicability of sparse BSS methods to real-world data, especially for largescale ones, largely depends on the design of reliable, efficient and versatile methods. Methods such as GMCA presenting these properties rely on heuristics that do not however guarantee any optimality of the estimate. On the contrary, the estimation using PALM algorithm to minimize a 1 regularized data fidelity term, which is now a standard approach to tackle generic matrix factorization problems due to mathematical guarantees, does not always exhibit such characteristics. To develop an optimization framework merging the best of the two worlds, we first investigated the behavior of PALM in the context of sparse BSS. These investigations show and explain the dramatic sensitivity of the estimate with respect to the regularization parameters (lack of efficiency and versatility) and also to the initialization (lack of reliability). To mitigate these limitations, we rationalize a hybrid approach that allows to circumvent the robustness issue of PALM with respect to the initialization and provides a proxy to automatically tune the regularization parameters for various data. Numerical experiments on both simulated and realistic data demonstrate the quality of the proposed approach with respect to the three desired characteristics. They also exhibits an improved separation accuracy with respect to state-of-the art methods. On the other hand, an extensive study of the limitations of the proposed strategy is performed, enlightening when to use it. As supplementary material, we also propose in Appendix E several exploratory alternative ways of how to perform a two-step strategy, which could pave the way for various enhancements.

(a) C d = 1 (b) C d = 2 (c) C d = 5 α (d) C d = 10 (e) C d = 20 (f) C d = 100 α (g) C d = 1 (h) C d = 2 (i) C d = 5 SNR (dB) α (j) C d = 10 SNR (dB) (k) C d = 20

Chapitre V

Tackling a high number of sources : blockGMCA

The goal of the previous chapter was to find both a reasonable initialization for a sparse BSS algorithm and a way to set the regularization parameters, which was a first mandatory step to hope to handle large-scale problems. In this chapter, we now enter the core subject of this thesis with a first large-scale issue of BSS : handling a large number of sources n. This problem is both difficult, since most methods fail when the number of sources typically exceeds a few tens, and of paramount importance in various applications such as spectroscopy, astronomy [Bobin et al. 2008] or biomedical imaging [Biswal & Ulmer 1999]. The proposed approach focuses on the optimization strategy, which has already been deemed as crucial in the previous chapters. More specifically, the proposed block-Generalized Morphological Component Analysis (bGMCA) algorithm builds upon block-coordinate descent with intermediate size blocks. Numerical experiments are provided that show the quality of the approach when the sources are numerous.

A Problem and outline

A.1 Problem : decreased performances with large numbers of sources

We illustrate in Fig. V.1 the performance deterioration of most BSS methods when the number of sources n becomes large. The evolution of the mixing matrix criterion as a function of the number of sources shows that most methods do not perform correctly in such a large-scale regime. In this case, the main source of deterioration is very likely related to the non-convex nature of BSS problem and a regularization issue. Indeed, for a fixed number of samples t, an increasing number of sources n will make these algorithms more prone to be trapped in spurious local minima, which tends to hinder the applicability of BSS on practical issues with a large n. Consequently, the optimization strategy has a huge impact on the separation performances.

A.2 Outline

The goal of this chapter is to introduce a novel algorithm dubbed bGMCA to specifically tackle sparse BSS problems when many sources need to be estimated. In addition to sparse modelling, this algorithm builds upon an efficient minimization scheme based on block-coordinate descent, as explained in Section B. In contrast to state-of-the art methods [Zibulevsky 2003,Bobin et al. 2015,Rapin et al. 2014[START_REF] Gillis | [END_REF], we show that block-based minimization with intermediate block sizes allows the bGMCA to dramatically enhance the separation performances for large n. This is demonstrated through comparisons with state-of-the art methods in Section C, which have been carried out on various simulation scenarios. The last part of the Chapter shows the flexibility of bGMCA, with an application to sparse and non-negative BSS in the context of spectroscopy.

B Proposed approach : use of intermediate block-sizes

B.1 State-of-art

As explained in Chapter III-C, to bypass both the coupling between A and S and the non-convexity of problem II.7, a common idea of several strategies (BCD [Tseng 2001], PALM [START_REF] Bolte | [END_REF]], ALS) is to benefit from the multi-convex structure of (II.7) by using blocks [Xu & Yin 2014] in which each sub-problem is convex. The minimization is then performed alternately with respect to one of the coordinate blocks while the other coordinates stay fixed, which entails solving a sequence of convex optimization problems. Most of the already existing methods can then be categorized in one of two families, depending on the block sizes :

-Hierarchical or deflation methods : these algorithms use a block of size 1. For instance, Hierarchical ALS (HALS) ( [START_REF] Gillis | [END_REF] and references therein, [Comon & Jutten 2010]) updates only one specific column of  and one specific row of Ŝ at each iteration. The main advantage of this family is that each subproblem is often much simpler as their minimizer generally admits a closed-form expression. Moreover, the matrices involved being small, the computation time is much lower. The drawback is however that the errors on some sources/mixing matrix columns propagate from one iteration to the other since they are updated independently.

-Full-size blocks : these algorithms use as blocks the whole matrices  and Ŝ

(the block size is thus equal to n). For instance, GMCA [Bobin et al. 2008] is part of this family. One problem compared to hierarchical or deflation methods is that the problem is more complex due to the simultaneous estimation of a high number of sources. Moreover, the computational cost increases quickly with the number of sources, since large-size problems need to be handled (e.g. in GMCA, the pseudo-inverse of the whole matrices  and Ŝ needs to be computed at each iteration).

The gist of the proposed bGMCA algorithm is to adopt an alternative approach that uses intermediate block sizes. The underlying intuition is that using blocks of intermediate size can be recast as relatively small-scale source separation problems, which are simpler to solve as testified by Fig. V.1. As a byproduct, these subproblems are also less costly to tackle. On the other hand, they are not small enough to incur dramatic error propagation.

B.2 Proposed approach

In the following, bGMCA minimizes the problem in eq. (II.7) with blocks, which are indexed by a set of indices I of size r, 1 r n. In practice, the minimization is performed at each iteration on submatrices of  (keeping only the columns indexed by I) and Ŝ (keeping only the rows indexed by I). We thereafter re-used the work described in the previous chapter and used a 2step minimization strategy using intermediate-size coordinate blocks. As before, it comprehends a GMCA as warm-up stage and PALM as refining stage, enabling the various benefits described in Chapter III. In the following, we describe both steps.

B.2.1 Warm-up stage

In the framework of the proposed bGMCA algorithm, the GMCA-based warmup stage uses blocks of size 1 r n and alternates between the update of some submatrices of  and Ŝ (these submatrices will be noted ÂI and ŜI ). The whole stage is summarized below : 0 -Initialize the algorithm with random Â(0) .

For each iteration (l) :

1 -A submatrix ŜI is now updated instead of the whole Ŝ. This is performed using a projected least square solution :

Ŝ(l) I = prox G(.) ( ÂI (l-1) † R I ) (V.1)
where : R I is the residual term defined by R I = X -ÂI C (l-1 ) Ŝ(l-1)

I C
(with I C the indices of the sources outside the block), which is the part of X to be explained by the sources in the current block I .

2 -The mixing sub-matrix ÂI is similarly updated with a fixed S :

ÂI (l) = prox J (.) (R I S (l) † I ) (V.2)
In this chapter, the penalizations J and G we will consider are the ones described in Section III-B.2.

B.2.2 Refinement stage

The PALM-based refinement stage using intermediate-size blocks reads as :

While the stopping criterion ∆ (l) has not reached the desired value, iterate over (l) :

1 -Update of a submatrix ŜI instead of the whole Ŝ :

Ŝ(l) I = prox γG(.) ÂI (l-1)T ÂI (l-1) 2  Ŝ (l-1) I - γ ÂI (l-1)T ÂI (l-1) 2 ÂI (l-1)T ( Â(l-1) Ŝ(l-1) -X)   (V.3) 2 -Update of a submatrix ÂI instead of the whole  : Â(l) I = prox δJ (.) Ŝ(l) I Ŝ(l) T I 2   ÂI (l-1) - δ Ŝ(l) I Ŝ(l) T I 2 ( Â(l-1) Ŝ(l) -X) Ŝ(l) T I   (V.4) 3 -Update stopping criterion : ∆ (l) = j∈[1,n] Âj (l) Âj (l-1) 1 n
Where the notations and the different constants are the same as in Chapter III.

B. Proposed approach : use of intermediate block-sizes 91

B.2.3 Block choice

Several strategies for selecting at each iteration new block indices I have been investigated :

-Sequential : at each iteration, r sources are selected sequentially in a cyclic way ;

-Random : at each iteration, r indices in [1, n] are randomly chosen following a uniform distribution and the corresponding sources updated ;

-Random sequential : this strategy combines the sequential and the random choices to ensure that all sources are updated an equal number of times.

In the experiments, random strategies tended to provide better results. Indeed, compared to a sequential choice, randomness is likely to make the algorithm more robust with respect to spurious local minima. Since the results between the random strategy and the random sequential one are similar, the first was eventually selected.

B.2.4 Convergence

The use of intermediate-size blocks do not alter the convergence guarantees of PALM algorithm with fixed thresholds. As such, the refinement stage converges to a stationary point of eq. (II.7), as long as the blocks are updated following an essentially cyclic rule [Chouzenoux et al. 2016] or even if they are chosen randomly and updated one by one [Patrascu & Necoara 2015].

B.2.5 Complexity

In this part, we focus only on the warm-up stage, which iterations are the most computationally expensive. Each iteration can then be decomposed into the following elementary steps : i) a residual term is computed with a complexity of O(mtr), where m is the number of observations, t the number of samples and r the block size ; ii) the pseudo-inverse is performed with the singular value decomposition of a r × r matrix, which yield an overall complexity of O(r 3 + r 2 m + m 2 r) ; iii) the thresholding-strategy first requires the evaluation of the threshold values, which has a complexity of rt ; iv) then the soft-thresholding step which has complexity O(rt) ; and v) updating A is finally performed using a conjugate gradient algorithm, whose complexity is known to depend on the number of non-zero entries in S and on the condition of this matrix C d (S). An upperbound for this complexity is thus O(rt C d (S)). The final estimate of the complexity of a single iteration is finally given by : r

[mt + rm + m 2 + r 2 + t C d (S)] (V.5)
With C d (S) the conditioning number of S. Thus, both the r factor and the behavior in r 3 show that small r values will lower the computational budget of each iteration. Since the algorithm is iterative, the final running time will however depend on both the complexity of each iteration and of the number of iterations. Intuitively, the required number of iterations should be inversely proportional to r, since only r sources are updated at each iteration, requiring n/r times the number of iterations needed by an algorithm using the full matrices. As will be emphasized later on, the number of required iterations will be smaller than expected, which on overall makes that the bGMCA algorithm enables a reduction of the computation time.

C Explaining the behavior of bGMCA : numerical experiments on simulated data

In this part, we present our results on simulated data. The goal is to show and to explain on simple data how bGMCA works.

C.1 Experimental protocol

The simulated data were generated in the following way :

1 -Source matrix S * : the sources are exactly sparse in the sample domain (that is, Φ S = I d -the results would however be identical for any source sparse in an orthogonal representation). Their coefficients are drawn randomly according to a Bernoulli-Gaussian distribution : among the t samples (t = 1 000), a proportion p (unless specified, p = 0.1) of the samples is non-zero, with an amplitude drawn according to a standard normal distribution.

2 -Mixing matrix A * : the mixing matrix is drawn randomly according to a standard normal distribution and modified to have unit columns and a given condition number C d (unless specified, C d = 1).

The number of observations m is taken equal to the number of sources : m = n. In this first simulation, no noise is added. The number of iterations for the warmup stage is 10 000. Here, the following penalizations G and J were used (for more mathematical details and the corresponding proximal operators, see Sec. III-B.2) :

-1 sparsity constraint in some transformed domain : The constraint G on S is a 1 -norm penalization.

-Oblique constraint : to avoid degenerated A and S matrices, the columns of A are constrained through J to lie onto the 2 hyper-sphere.

To measure the accuracy of the separation, we again followed the definition in [Bobin et al. 2015] to use the global criterion C A on A (cf. Appendix B). The data matrices being drawn randomly, each experiment was performed several times (typically 25 times) and the median of C A over the experiments will be displayed.

C.2 Modeling block minimization

In this section, a simple model is introduced to describe the behavior of the warm-up stage of the bGMCA algorithm. As described in section B.2, updating a given block is performed at each iteration from the residual R I = X -ÂI C ŜI C (to C. Explaining the behavior of bGMCA : numerical experiments on simulated data 93 lighten the notations, we dropped the iteration number l -1). If the estimation were perfect, the residual would be equal to the part of the data explained by the true sources in the current block indexed by I, which would read : R I = A I * S * I . It is nevertheless mandatory to take into account the noise N, as well as a variety of flaws in the estimation by adding a term E to model the estimation error. This entails :

R I = X -ÂI C ŜI C = A I * S * I + E + N (V.6)
A way to further describe the structure of E is to decompose the estimated ŜI matrix in the true matrix plus an error : ŜI = S * I + s I and ŜI C = S * I C + s I C , where s is the error on S * . Assuming that the errors are small and neglecting the second-order terms, the residual R I can now be written as :

R I = X -ÂI C ŜI C = A I * S * I + A I C * S * I C -ÂI C Ŝ * I C -ÂI C s I C + N (V.7)
This implies that :

E = (A I C * -ÂI C )S * I C -ÂI C s I C (V.8)
Equation (V.8) highlights two terms. The first term can be qualified as interferences in that it comes from a leakage of the true sources that are outside the currently updated block. This term vanishes when ÂI C is perfectly estimated. The second term corresponds to interferences as well as artefacts. It originates indeed from the error on the sources outside the block I. The artefacts comprehend in particular the errors on the sources induced by the soft thresholding corresponding to the 1 -norm. Equation (V.8) also allows us to understand how the choice of a given block size r n will impact the separation process :

-Updating small-size blocks can be recast as a small-size source separation problem where the actual number of sources is equal to r. As testified by Fig.

V.1, updating small-size block problems should be easier to tackle.

-Small-size blocks should also yield larger errors E. It is intuitively due to the fact that many potentially badly estimated sources in I C are used for the estimation of A I * and S * I through the residual, deteriorating this estimation. It can be explained in more details using equation (V.

C.3 Experiment

In this section, we investigate the behavior of the proposed block-based GMCA algorithm with respect to various parameters such as the block size r, the number of sources n, the conditioning of the mixing matrix C d and the sparsity level of the sources p.

C.3.1 Study of the impact of r and n

In this subsection, bGMCA is evaluated for different numbers of sources n = 20, 50, 100. Each time the block sizes vary in the range 1 ≤ r ≤ n. In this experiment and to complete the description of section C.1, the parameters for the matrix generation were : p = 0.1, t = 1 000, C d = 1, m = n, with a Bernoulli-Gaussian distribution for the sources. These results are displayed in Fig. V.2a. Interestingly, three different regimes characterize the behavior of the bGMCA algorithm :

-For intermediate and relatively large block sizes (typically r > 5 and r < n-5) :

we first observe that after an initial deterioration around r = 5 , the separation quality does not vary significantly for increasing block sizes. A degradation of several dB can then be observed for r close to n. In all this part of the curve, the error term E is composed of residuals of sparse sources, and thus E will be rather sparse when the block size is large. Based on the MAD, the thresholds are set according to dense and not to sparse noise. Consequently the automatic thresholding strategy of the bGMCA algorithm will not be sensitive to the estimation errors.

-A very prominent peak can be observed when the block size is of the order of 3. Interestingly, the maximum yields a mixing matrix criterion of about 10 -16 , which means that perfect separation is reached up to numerical errors. This value of 160 dB is at least 80 dB larger than in the standard case r = n, for which the values for the different n are all below 80 dB. In this regime, error propagation is composed of the mixture of a larger number of sparse sources, which eventually entails a densely distributed contribution that can be measured by the MAD-based thresholding procedure. Therefore, the threshold used to estimate the sources is able to filter out both the noise and the estimation errors. Moreover, r = 5 is quite small compared to n. Following the modeling introduced in section C.2, small block sizes can be recast as a sequence of low-dimensional blind source separation problems, which are simpler to solve.

-For small block sizes (typically r < 4), the separation quality is deteriorated when the block size decreases, especially for large n values. In this regime, the level of estimation error E becomes large, which entails large values for the thresholds Λ. Consequently, the bias induced by the soft-thresholding operator increases, which eventually hampers the performance quality. Furthermore, for a fixed block size r, E increases with the number of sources n, making this phenomenon more pronounced for higher n values.

C.3.2 Condition number of the mixing matrix

In this section, we investigate the role played by the conditioning of the mixing matrix on the performances of the bGMCA algorithm. Gaussian distribution and p = 0.1, t = 1000. One can observe that when C d increases, the peak present for r close to 5 tends to be flattened, which is probably due to higher projection errors. At some iteration l, the sources are estimated by projecting X-ÂI C ŜI C onto the subspace spanned by ÂI . In the orthogonal case, the projection error is low since ÂI C and ÂI are close to orthogonality at the solution.

However, this error increases with the condition number C d .

C.3.3 Sparsity level p

In this section, the impact of the sparsity level of the sources is investigated. The sources are still following a Bernoulli-Gaussian distribution. The parameters are : n = 50, t = 1 000, C d = 1. As featured in Figure V.3, the separation performances at the maximum value decrease slightly with larger p, while a slow shift of the transition between the small/large block size regimes towards larger block sizes operates. Furthermore, the results tend to deteriorate quickly for small block sizes (r < 4). Indeed, owing to the model of subsection C.2, the contribution of S * I C and s I C in the error term (V.8) increases with p, this effect being even more important for small r (which could also explain the shift of the peak for p = 0.3, by a deterioration of the results at its beginning, r = 3). When p increases, the sources in ŜI become denser. Instead of being mainly sensitive to the noise and E, the MAD-based thresholding tends to be perturbated by ŜI , resulting in more artefacts, which eventually hampers the separation performances. This effect increases when the sparsity level of the sources decreases.

C.3.4 Number of iterations and computation time

We have already seen in Section B.2.5 that the bGMCA algorithm enables a gain in terms of computational complexity of each iteration (roughly speaking, this gain is almost linear if t is much larger than m and r). We here further empirically assess the actual number of iterations required by the warm-up stage to yield a good initialization. To this end, the following experiment has been conducted : 1. First, the algorithm is launched with a large number of iterations (e.g. 10 000)

to give good  and Ŝ matrices. The corresponding value of C A is saved and called C * A . 2. Using the same initial conditions, the warm-up stage is re-launched and stops when the mixing matrix criterion reaches 1.05 × C * A (i.e. 5% of the "optimal" solution for a given setting).

The number of iterations needed to reach the 5% accuracy is reported in Fig. V.4. Intuitively, one would expect that when the block size decreases, the required number of iterations should increase by about n/r to keep the number of updates per source constant. This trend is displayed with the straight green curve of Fig. V.4. Interestingly,Fig. V.4 shows that the actual number of iterations to reach the 5% accuracy criterion almost does not vary with r. Consequently, on top of leading to computationally cheaper iterations, using small block sizes almost requires the same number of iterations for the warm-up stage to give a good initialization. Therefore, the use of blocks allows a huge decrease of the computational cost of the warm-up stage and thus of sparse BSS.

D Validation of the approach on realistic sources

D.1 Context

The goal of this part is to evaluate the behavior of bGMCA and show its efficiency in a more realistic setting. Our data come from a simulated LC -1 H NMR (Liquid Chromatography -1 H Nuclear Magnetic Resonance) experiment. The objective of such a experiment is to identify each of the chemicals compounds present in a fluid, as well as their concentrations. As explained in Chapter II-A.1.2 (the principles of LC -1 H NMR and LC -MS experiments are the same), the LC -1 H NMR experiment enables a first physical imperfect separation during which the fluid goes through a chromatography column and its chemicals are separated according to their speeds (which themselves depend on their physical properties). Then, the spectrum of the output of the column is measured at a given time frequency. These measurements of the spectra at different times can be used to feed a bGMCA algorithm to refine the imperfect physical separation. The fluids on which we worked could for instance correspond to drinks. The goal of bGMCA is then to identify the spectra of each compound (e.g. caffeine, sucrose, menthone...) and the mixing coefficients (which are proportional to their concentrations) from the LC -1 H NMR data. BSS has already been successfully applied [Toumi et al. 2013] to similar problems but generally with lower number of sources n. The sources (n = 40 sources with each t = 10 000 samples) are composed of elementary sparse non-negative theoretical spectra of chemical compounds taken from the SDBS database 1 , which are further convolved with a Laplacian having a width of 3 samples to simulate a given spectral resolution. Therefore, each convolved source becomes an approximately sparse non-negative row of S * (cf.

D.2 Experiments

There are two main differences with the previous experiments of section C : i) the sources are sparse in the undecimated wavelet domain Φ S , which is chosen as the starlet transform [START_REF] Starck | [END_REF] in the following, and ii) the non-negativity of S and A is enforced. These results show that non-negativity yields a huge improvement for all block sizes r, which is expected since the problem is more constrained. This is probably due to the fact that all the small negative coefficients are set to 0, thus artificially allowing lower thresholds and therefore less artefacts. This is especially advantageous in the present context with very low noise2 (the Signal to Noise Ratio -SNR -has a value of 120 dB) where the thresholds do not need to be high to remove noise. Furthermore, the separation quality tends to be constant for r ≥ 10. In this particular setting, non-negativity helps curing the failure of sparse BSS when large blocks are used. However, using smaller block sizes still allows reducing the computation cost while preserving the separation quality. The bGMCA with non-negativity also compares favorably with respect to other tested standard BSS methods, yielding better results for all values of r. A single original source is displayed in the right panel of Fig. V.6 after its convolution with a Laplacian. Its estimation using bGMCA with a non-negativity constraint is plotted in dashed line on the same graph, showing the high separation quality because of the nearly perfect overlap between the two curves. Both sources are drawn in the direct domain.

The robustness of the bGMCA algorithm with respect to additive Gaussian noise has further been tested. Fig. V.7 reports the evolution of the mixing matrix criterion for varying values of the signal-to-noise ratio. It can be observed that bGMCA yields the best performances for all values of SNR. Although it seems to particularly benefit from high SNR compared to HALS and EFICA, it still yields better results than the other algorithms for low SNR despite the small block size used (r = 10), which could have been particularly prone to error propagations.

SNR (logarithmic scale)

C A (dB)

Conclusion

While being central in numerous applications, tackling sparse BSS problems when the number of sources is large is highly challenging. In this chapter, we described the block-GMCA algorithm, which is specifically tailored to solve sparse BSS in such a large-scale regime. In contrast to other state-of-the-art algorithms, bGMCA builds upon block-coordinate optimization with intermediate-size blocks. The mechanisms enabling bGMCA to have improved performances compared to its full-size block counterparts are explained both using exactly sparse simulated data and a mathematical modeling. While on such exactly sparse data bGMCA can lead to numerically perfect separations, comparisons have also been carried on simulated spectroscopic data, which demonstrates the reliability of the proposed algorithm in a realistic setting and its superior performances for high SNR. All the numerical comparisons conducted show that bGMCA performs at least as well as standard sparse BSS on mixtures of a high number of sources and most of the experiments even show dramatically enhanced separation performances. As a byproduct, the proposed block-based strategy yields a significant decrease of the computational cost of the separation process.

Chapitre VI

Tackling large-scale datasets : mini-batch optimization with aggregation on Riemannian manifold

The above chapter proposed a way to handle problems involving a high number of sources n, that is matrix factorization problems with a large inner dimension. In this chapter, we introduce a scalable sparse BSS algorithm enabling to handle large-scale datasets X, which is of uttermost importance due to current ever growing data-sizes. The proposed distributedGMCA (dGMCA) combines a robust projected alternating least-squares method with mini-batches, which enables to both handle large-scale datasets and to benefit from the efficiency, reliability and versatility of pALS (see Chapter IV). The originality lies in the use of a manifold-based aggregation of the different estimates of the mixing matrices. This approach is showed to maintain high performances compared to full-batch methods, which are in contrast not distributed. Remarkably, dGMCA can further outperform such algorithms when the sources have highly sparse distributions. Numerical experiments are carried out on synthetic data as well as realistic simulations of spectroscopic data.

This work has been started by the internship of Tobias Liaudat, which was proposed during the present PhD.

A Introduction

As explained in Chapter II, the challenges implied by huge datasets X are twofold : time computational issues (the dataset sizes make that it is currently impossible to perform BSS in a decent amount of time) and memory limitation issues (datasets are even too large to fit into memory). To give an order of the magnitude at stake, one could expect from astronomical devices such as the SKA data of up to m 10 4 observations and t 10 9 samples. In this context, the GMCA algorithm [Bobin et al. 2007, Bobin et al. 2015], which was shown in Chapter IV to perform well for handling small-scale to middle [Bobin with aggregation on Riemannian manifold et al. 2013] size datasets X, is not anymore usable. Indeed, it needs at each iteration both an inversion of the factors A and S and multiplications with X, which is particularly costly. This is furthermore an issue for the two-step approach of Chapter IV, that at least required a few iterations of GMCA as a warm-up stage for PALM. In this chapter, we will focus on the largest dimension of X : we will propose a new sparse BSS method enabling to cope with datasets X comprehending a large number of samples t. Before detailing the method, we will present the context by reviewing other related works. Then, we explain the challenges and our contribution.

A.1 Sparse matrix factorization for large-scale datasets X

The sparse BSS problem of Eq. (II.8) can be seen as a generic sparse matrix factorization problem, for which some works have been dedicated to large-scale datasets. A classical idea to tackle such an issue is to use only a submatrix (a set of columns) of X at each iteration, that is to use mini-batches. In contrast to successful small scale sparse BSS algorithms, most of the corresponding works do not use pALS. Among them, one can distinguish the ones that :

-builds on stochastic gradient descent (SGD - [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF]]) : for instance, [Davis et al. 2016] extended the use of PALM to mini-batches, making it possible to tackle large datasets. Such a (potentially asynchronous) approach has also been extended to the case of hyperspectral imaging [Thouvenin et al. 2018]. In this work, the authors use the framework of [Cannelli et al. 2016], arguing a higher flexibility than in [Davis et al. 2016] in which the asynchronicity has a high impact on the allowable step sizes, counter-balancing its positive effects ;

-builds upon stochastic approximations. In the context of dictionary learning, [Mairal et al. 2009] proposed an online algorithm in which the dictionary is computed only 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 t b values. More recently, this algorithm has been extended in [Mensch et al. 2018] to tackle datasets that, in our context, would be huge both in the number of samples t and number of observations m. An extension has also been envisioned in the context of hyperspectral imaging with spectral variability [START_REF] Thouvenin | Online unmixing of multitemporal hyperspectral images accounting for spectral variability[END_REF].

A.2 Challenges and contributions

A.2.1 Challenges

While the approaches of the previous subsection using mini-batches might sound appealing for large-scale sparse BSS, it must be highlighted that they have been elaborated for different kind of problems. In these, they generally do not aim at retrieving physical factors A * and S * , which might not even exist. As such, algorithms based on stochastic approximations as in [Mairal et al. 2009,Mensch et al. 2018] were initially targeting the dictionary learning problem and thus focus only on finding factors yielding good results for a given application (e.g. denoising...). Therefore, they generally do not provide accurate results in the context of sparse BSS (this claim will be further backed by the comparisons performed throughout this whole work). On the other hand, it has been shown in Chapter IV that gradient descent (GD) methods usually suffer from a low reliability in the sparse BSS context. Indeed, and in contrast to the pALS scheme of [Bobin et al. 2007], the use of GD methods makes automatic hyper-parameter choice much more difficult. Therefore, the mini-batches SGD methods are not expected to work well either (at least currently without requiring a warm-up stage based on pALS). This leaves sparse BSS with very few satisfying options, as :

-Large-scale optimization algorithms using mini-batches yield bad physical factors  and Ŝ in practice ;

-Usual sparse BSS algorithms working well on small-scale problems are not scalable.

Consequently, the solution we propose for large-scale BSS is to introduce minibatches in the GMCA algorithm of [Bobin et al. 2007], which would enable to benefit from the best of the two worlds : re-using the mini-batch approach enabling scalable algorithms ; benefiting from the automatic parameter choice yielded by the pALS scheme and its reliability.

A.2.2 Other related works

The ALS algorithm has already been studied in the large-scale setting, but in the different application of recommendation systems. In [START_REF] Zhou | [END_REF], the authors proposed a parallel ALS with weighted-λ-regularization to solve a low-rank matrix factorization problem. To do that, they split the X matrix, the factors A and S accordingly, and each node processes its own submatrices, in an approach reminiscent to mini-batches. The method was extended by [Teflioudi et al. 2012] to a distributed (shared nothing) setting. The works of [START_REF] Hastie | [END_REF][START_REF] Christian | Parallelization of the alternating-least-squares algorithm with weighted regularization for efficient gpu execution in recommender systems[END_REF]] also tackle a similar issue. However, several differences between such a problem and ours must be highlighted as : i) The application is different from ours : in particular, as X is sparse (due to the high number of missing entries), it is much easier to work on smaller submatrices and the mini-batches are much more naturally chosen ; ii) As a consequence, no aggregation (see below) is used since during each iteration each matrix  and Ŝ (or more specifically, a submatrix of  and Ŝ) is updated once ; iii) The cost function is different from ours, as there is no 1 sparsity promoting term. As such, while the main idea of using mini-batches within an ALS scheme will be re-used in this Chapter, the similarities with these works are limited. with aggregation on Riemannian manifold

A.2.3 Contributions

In this Chapter, we propose to extend sparse BSS to the large-scale setting, enabling to tackle datasets X with a huge number of columns t. Our approach consists in introducing a mini-batch (stochastic) version of [Bobin et al. 2007] coined distri-butedGMCA (dGMCA). The algorithm enables to tackle datasets that would not even fit into the memory of one computer. The approach is furthermore robustified through the introduction of an aggregation step based on a robust weighted mean on a Riemannian manifold, enabling to aggregate several estimators of the mixing matrix A. Beyond enabling scalability, the introduction of such a step is empirically shown to improve in some experiments the results of dGMCA over the full batch version of the algorithm. Lastly, realistic experiments are carried out to demonstrate the quality of the approach, and relationship to SGD methods in machine learning are highlighted to explain the results.

B Distributed sparse Alternating Least-Squares

B.1 Distributing the GMCA algorithm

B.1.1 Naive approach

A simple approach to deal with large-scale data would be to split the large dataset X into B submatrices. As such, the t columns would be split into B disjoint submatrices, each of them having t b = t/B columns (for the moment, t is assumed to be a multiple of B). Each of the corresponding submatrix is denoted as

X J b , b ∈ [1, B]
and the corresponding indices of the columns as J b , b ∈ [0, B] (#J b = t b ). A naive approach would then be to work independently on each submatrix X J b , b ∈ [1, B]. As such, instead of looking for an (approximate) minimizer of Eq. (II.8) using GMCA, we would rather tackle small-scale subproblems :

argmin A[J b ]∈R m×n ,S J b ∈R n×t b 1 2 X J b -A[J b ]S J b 2 F + R S J b S J b 1 +ι {∀i∈[1,n], A[J b ] i 2 2 =1} (A[J b ]), (VI.1)
where, for the sake of simplicity, the Φ S has been omitted and taken equal to identity. The final estimate of S * is straightforwardly obtained through the concatenation of the columns of the different final estimates

ŜJ b , b ∈ [1, B].
The question is however more intricate concerning the mixing matrix, as each submatrix X J b yields a different estimate A[J b ] of the same full matrix A * . At this point, we should make clear that the A[J b ] are different despite the assumed stationarity of the signals, as they rely on a specific small-case realization X J b . As such, some will be better estimates of the true underlying A * than others. To get the final estimate Â, we thus need to resort to an aggregation step that will merge the information yielded by the various A[J b ] to yield a final  :

 = AGGREGATE(A[J 1 ], A[J 2 ], ..., A[J B ]), (VI.2)
where the AGGREGATE function is used as a generic term for the aggregation. As a simple example of such a function, one could for instance use the coefficient-wise Euclidean mean of the different

A[J b ], b ∈ [1, B]
. This option is however simplistic and can be enhanced (cf. Sec. B.2). The naive version of dGMCA is summarized in Algorithm 5.

Algorithm 5 Naive dGMCA

1: procedure NAIVE dGMCA(X, Â(0) ) 2: Choose J 1 , J 2 , • • • , J B as a partition of [1, t] 3: for b = 1, • • • , B do 4: A[J b ], ŜJ b = GMCA(X J b , Â(0) ) 5:
end for

6: Â = AGGREGATE(A[J 1 ], A[J 2 ], ..., A[J B ]) 7:
return Â, Ŝ 8: end procedure Nevertheless, one of the main flaws of the naive dGMCA is that the information yielded by a given submatrix X J b 1 is not shared at all to enhance the mixing matrix A[J b ] found from another submatrix X J b , b = b 1 . Said differently, processing each submatrix X J b fully independently makes that much less information than using the whole matrix X is used to estimate each of the A[J b ]. It could be argued that this is not a real issue, as we are not directly interested by each of the A

[J b ], b ∈ [1, B],
but rather by the final  obtained after the aggregation step, which goal is precisely to merge the information between the subproblems. Unfortunately, this is not fully accurate and practical experiments further tend to hinder such arguments. Indeed, performing independent GMCA makes that none of the estimated A[J b ] is descent (as only too small submatrices of the original data X are used). As such, the aggregation step -said differently, information sharing among subproblems -is applied "too late" on only badly estimated A[J b ] : since  is estimated from mainly bad A[J b ], it is also bad.

B.1.2 dGMCA algorithm

To alleviate the previous difficulty induced by working independently on subproblems of the form Eq. (VI.1), we propose to enable the estimation of each A[J b 1 ] to share information with the other A[J b ], b = b 1 during the iterations of GMCA. To do that, we propose to apply the aggregation step during each iteration of GMCA (following the machine learning terminology, we can now speak of epochs of GMCA), in contrast to the naive approach in which it was performed only after the end of GMCA. Using such an approach, the indices J b may now change during GMCA iterations, as well as the corresponding submatrices X J b . As such, we will rather write for each epoch J b (l). By analogy with the works of [Mairal et al. 2009, Davis et al. 2016, with aggregation on Riemannian manifold Mensch et al. 2018], we will further call X J b (l) a mini-batch of X 1 . A last remaining question before detailing the algorithm is how to choose the regularization parameters R S within dGMCA, to keep the advantage of the adaptive parameter tuning of GMCA 2 . Indeed, we will here aim at extending the most recent parameter choice of GMCA 3 , based on using an increasing percentile of the whole distribution of the currently estimated Ŝ(l) , which is difficult to distribute for large t values. As such, we propose to rather use an exponential decay of the parameters of the form :

R i(l) i S i = κσ i + Ŝ(l) i ∞ -κσ i exp(-lα i ), (VI.3) where R i(l) i S i
is the estimated regularization parameter at epoch l for source S i , σ i is an estimation of the noise back-projected on source S i , κ is chosen according to the fixed point Gaussian noise removal argument (here κ = 3), Ŝ(l)

i ∞ is the maximum absolute value of Ŝ(l) i and α i is a parameter controlling the exponential decay decrease 4 . An interesting property of the parameter choice of (VI.3) is that it is highly distributable :

Ŝ(l) i ∞
is the maximum over each mini-batch

( Ŝ(l) i ∞ = max b∈[1,B] ŜJ b (l) i ∞
), and σ i can be chosen as the median of its estimations over the mini-batches :

σ i = median b∈[1,b] σ i [J b (l)]. Each estimation σ i [J b (l)]
is in turn performed according to the usual method used in GMCA, based on the MAD operator (cf. Chapter II [Bobin et al. 2007]).

B.1.3 Summary of the algorithm

B.2 Manifold-based mixing matrix aggregation

We now detail the aggregation step used to build the estimate Â(l) from the various

A (l) [J b ], b ∈ [1, B].
To lighten the notations, we will drop in this subsection the l index, but the aggregation is of course performed at each iteration. We propose to define  as the barycenter of the different estimates according to some φ as follows :

1. Note that in these works, the aggregation step enabling to share the information between mini-batches is however implicit and made simpler due to the fact that the algorithms use only local updates such as gradient steps instead of least-square solutions.

2. While in the naive dGMCA approach the answer to such a question was straightforward, it is important to emphasize that choosing the regularization parameters using small X J b could lead to estimate them badly.

3. The reader might wonder why we are not implementing here the regularization parameter choice used in Chapter IV, which was fully based on the MAD. This is due to the fact that while this would is easier to distribute, the use of the percentile decrease generally provides better results in GMCA. In contrast, the experiments we perform in the Appendix E show that the percentile does not seem to work well within PALM, thus our choice to use the MAD.

4. The αi parameter intrinsically depends on the sparsity level of Si. As such, it can be estimated by fitting a generalized Gaussian distribution to the current estimation Ŝ(l) i and using a maximum likelihood estimator. However, we experimentally found out that the final result of dGMCA is quite robust to the estimation of αi, enabling to rather use a value fixed beforehand.

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

A[J b (l)] = Π . 2 =1 (X J b (l) ŜJ b (l) † ) 7:
end for 8:

Â(l) = AGGREGATE(A[J 1 (l)], A[J 2 (l)], ..., A[J B (l)]) 9:
end for 10:

return Â(L) , Ŝ(L) 11: end procedure

 = argmin A∈R m×n B b=1 ω b φ (A, A[J b ]) , (VI.4)
where the barycentric weights are positive and sum to one : ∀b ∈ [1, b]; ω b ≥ 0 and B b=1 ω b = 1. A straightforward example amounts to choose φ as the standard Euclidean distance applied on each column :

φ (A, A[J b ]) = n j=1 A j -A j [J b ] 2
2 . This choice will eventually define the aggregated estimator as a weighted sum of the different minibatch estimators :

 = B b=1 ω b A[J b ].
(VI.5) However, in the context of BSS, the mixing matrix is assumed to belong to the Oblique ensemble. In the limit of small angular distances between the different estimators, the Euclidean metric is likely an aggregated estimate that fulfills the Oblique constraint, which can be further constrained by projecting  onto O b . In the general case, this is unlikely to hold true, especially when considering small size mini-batches, that can lead to more larger angular deviations. This implies that the Oblique constraint has to be preserved in the aggregation step.

Fréchet mean on the hypersphere

Recall that the Oblique constraint implies that each column of the mixing matrix A belongs to the m-dimensional hypersphere S m or m-sphere, which is a Riemannian manifold. We assume the reader to be familiar with some basic notions about optimization on Riemannian manifolds. If it is not the case, we refer to [START_REF] Absil | [END_REF] or to the Appendix F for a summary of useful elements for this work.

A natural approach to take into account the Oblique constraint consists in building an aggregated estimator by defining each of the estimate  column as being with aggregation on Riemannian manifold the Fréchet mean of the corresponding columns in the estimators A[J b ] :

∀j ∈ [1, n], Âj = argmin a∈R m b∈[1,B] ω b d β a, A[J b ] j , (VI.6)
where 1 ≤ β < +∞. The case β = 2 corresponds to the 2 norm along the geodesics of the manifold S m . Following [Afsari 2011], it is possible to find a local critical point of this problem using an iterative gradient descent algorithm on the m-sphere presented in Algorithm 7. The parameter ρ is the step size 5 . The term ∇J β ( Âj (l f ) ), where l f is the iteration number in the Fréchet mean algorithm, is the gradient of the mean cost function

J β ( Âj (l f ) ) = B b=1 ω b d β Âj (l f ) , A j [J b ]
, which takes the following expression for β ≥ 1 :

∇J β ( Âj (l f ) ) = - B b=1 ω b d β-2 Âj (l f ) , A j [J b ] log Âj (l f ) A j [J b ] . (VI.7) Algorithm 7 Fréchet mean 1: procedure Fréchet mean on the Oblique ensemble O b (m) 2:
Correct for permutations 3:

for j = 1, • • • , n do
Loop over all the sources 4:

while Convergence is not reached do 5:

∇J β ( Âj (l f ) ) = -B b=1 ω b log Âj (l f ) A j [J b ] Gradient 6: Âj (l f +1) = exp Âj (l f ) -ρ∇J β ( Âj (l f )
))

7:

l f ← l f + 1 8:
end while 9:

end for 10:

Return Â(L f ) 11: end procedure Robust Fréchet mean on the hypersphere The use of small mini-batches makes the separation process more sensitive to several factors. In practice, this will tend to generate outliers in the estimated A[J b ] for some mini-batches, which is discussed in more details in Section B.3 and supported by several numerical experiments in Section C. Unfortunately, the Fréchet mean is not robust to such outliers [Fletcher et al. 2008, Arnaudon et al. 2013]. In this case, a natural choice would be to choose β = 1, which corresponds to the usual 1 norm. However, the 1 norm is not differentiable about 0 and the gradient of the cost function depends on the inverse of the distance d 1 as highlighted in Equation (VI.7), which makes it quite unstable in practice. To alleviate this issue, 5. It can possibly vary during the optimization process. It will be kept fixed in this Chapter.

we propose to build a differentiable approximation of d 1 based on Nesterov's smoothing technique [Nesterov 2005]. Such a smooth approximation of d 1 can be built as follows [Becker et al. 2011] for a, b ∈ R m :

d 1 ν (a, b) = argmax u ∞≤1 < a -b, u > - ν 2 u 2 2 . (VI.8)
This approximated distance is differentiable and its gradient is ν-Lipschitz. This entails that the gradient of the cost function takes the form :

∇J β ( Âj (l f ) ) = - B b=1 ω p ∇d 1 ν Âj (l f ) , A j [J b ] log Âj (l f ) A j [J b ] , (VI.9)
where the gradient of d 1 ν is given by [Becker et al. 2011] for all i ∈ [1, m] : Correct for permutations 3:

∇d 1 ν ( Âj (l f ) ) i =    ν -1 Âj (l f ) i , if | Âj (l f ) i | < ν, sign( Âj (l f ) i ),
for j = 1, • • • , n do
Loop over all the sources 4:

while Convergence is not reached do 5:

∇J β ( Âj (l f ) ) = -b∈[1,B] ω j b ∇d 1 ν Âj (l f ) , A j [J b ] log Âj (l f ) A j [J b ] 6: Âj (l f +1) = exp Âj (l f ) -ρ∇J β ( Âj (l f ) ) 7: l f ← l f + 1 8:
end while 9:

end for 10:

Return Â(L f ) 11: end procedure

From Fréchet mean to barycenter

In practical settings, the data are assumed to be contaminated with Gaussian noise. As pointed out in the previous paragraph, the use of small size mini-batches might lead to very diverse estimates of the mixing matrix. These estimates could yield sources with largely different signal-to-noise ratio. A first solution could be to choose weights {ω b } b∈[1,B] that penalize mini-batches with smaller SNR. However this would not take into account the actual noise that would contaminate the estimated sources. Therefore, we propose to choose these weights as being a function of the SNR of the estimated sources. Assuming that the data noise covariance matrix with aggregation on Riemannian manifold is denoted Σ N , the weights are defined as :

∀b ∈ [1, B]; ω j b = W[J b ] j Σ N W[J b ] j T -1 B p=1 W[J b ] j Σ N W[J b ] j T -1 (VI.11)
where

W[J b ] j = [A[J b ] + ] j .
In the numerical experiments, this weighting strategy will be applied to both the Fréchet mean and its robust version.

B.3 Another point of view about dGMCA : connections with stochastic gradient descent

To give a better insight of dGMCA principle, we aim in this subsection to highlight some connections with works on SGD.

From GD to ALS

Beyond matrix factorization problems, GD is a very popular approach in machine learning. In the scope of BSS and forgetting the oblique constraint, estimating A * using GD would yield the following update rule at epoch l for estimate Ŝ of the sources :

Â(l+1) = Â(l) + η (l) ∆ (VI.12) = Â(l) + η (l) X -Â(l) Ŝ ŜT (VI.13)
where ∆ is the gradient of the data-fidelity term of Eq. (II.8) with respect to Ŝ and η (l) the learning rate (i.e. the gradient step size when using the machine learning terminology). On the other hand, it is possible to incorporate second order information, which amounts to write a Newton update as follows :

Â(l+1) = Â(l) + η (l) X -Â(l) Ŝ ŜT H -1 (VI.14) = Â(l) + η (l) X -Â(l) Ŝ ŜT ŜŜ T -1 (VI.15)
where

H = ŜŜ T -1
is the Hessian of the data fidelity term with respect to A. Assuming the sources to be decorrelated yields a diagonal Hessian, which is equivalent to a single iteration of the standard GD algorithm up to a scaling factor and with

η (l) = 1. More generally, fixing η (l) = 1 entails : Â(l+1) = X ŜT ŜŜ T -1
. In this setting, ALS can be regarded as a special type of GD when curvature or secondorder information is used.

From SGD to mini-batch ALS However, when it turns to mini-batch optimization the connections become slightly delicate : whether in Newton descent algorithm or in ALS, the Hessian matrix ŜJ b ŜJ b T becomes mini-batch dependent. Small size mini-batches entail extra errors on the estimated Hessian, which eventually leads to more stochasticity on the estimated mixing matrices. Therefore, ALS leads to more stochasticity on the estimates than standard SGD, which can be detrimental to the optimization procedure when the mini-batch size becomes very small 6 .

From averaged SGD to aggregated stochastic ALS

The parallel between mini-batch ALS and SGD can be prolonged to give an insight concerning the aggregation step. More precisely, it is well known that the stochasticity of SGD can lead to an optimization path with more fluctuations and decreased rates of convergence in comparison to full-batch GD updates. To bypass this issue, a classical modification of SGD is to average the iterates of the estimates (Polyak-Rupper averaging [Ruppert 1988, Polyak & Juditsky 1992]), improving the convergence rates and reducing the impact of noise. In such an averaging process, the estimate of a variable Â(l) at iteration l is chosen as the mean over the previous iterations : Â(l) = 1 l l i=1 Â(i) , which resembles the Euclidean aggregation step discussed in Section B.2. Recent works [Tripuraneni et al. 2018] have extended the Polyak-Rupper averaging to a Riemannian setting, creating links with the aggregation step we use in dGMCA : instead of using a mere average of the previous iterates, the averaging is performed on a manifold. This parallel must however be tempered as : i) we apply a weighted Polyak-Rupper averaging : in particular, the weights are zero for the estimates found before the current epoch l 7 ; ii) we do not aggregate after each mini-batch, but rather once per epoch, which helps distributing the algorithm ; iii) [Tripuraneni et al. 2018] explores SGD, and we deal with pALS, implying the differences evoked above. However, we highlight that as pALS might increase the discrepancy between the different estimations compared to GD, such an aggregation might be particularly relevant to smooth the optimization path. As a side remark, while more different from the dGMCA aggregation, we would like to point out the works of [START_REF] Sato | [END_REF],Zhang et al. 2016], which are also variance reduced Riemannian optimization methods.

C Numerical experiments

In the following experiments, the performances of the dGMCA algorithm are evaluated in various experimental settings, on both simulated and realistic sources.

6. On the other hand, we shall see later that such stochasticity can also have benefits to some extent.

7. Note that using specific weights has been proved in some different settings to lead to interesting results through creating an implicit regularization [Neu & Rosasco 2018]. with aggregation on Riemannian manifold

C.1 Experiments on simulated data

Comparison set-up

In this subsection, we first make use of synthetic random data, which allows to perform Monte-Carlo simulations to assess the robustness of the different methods and to study the performances of dGMCA when varying the experimental parameters. More precisely, we will look at the influence of different numbers of sources n and observations m, as well as the one of mixing matrices A * with different condition numbers and sources with various sparsity levels p. To that end, the data are synthesized as follows :

-The sources S * j , j ∈ [1, n] have entries which are distributed independently and identically according to a Generalized Gaussian distribution with parameter 0 < α ≤ 1.

-The mixing matrix A * is picked at random from a Gaussian distribution, and further processed to have columns with unit 2 norm and a pre-defined condition number.

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

Beyond the proposed two versions of the dGMCA algorithm, namely the one using as aggregation the Fréchet mean and the one using the robust Fréchet mean, comparisons will be carried out with :

-GMCA : This algorithm is used as a baseline to compare its parallelized dGMCA counterpart.

-Online dictionary learning : see Section A.1 and [Mairal et al. 2009]. This algorithm is a classical one for solving large-scale sparse matrix factorization problems.8

To assess the separation quality, we will use the mixing matrix criterion C A .

Studying the impact of the number of observations m

In this paragraph, we evaluate the performances of the dGMCA with respect to the number of observations m. The role played by this parameter is twofold. First, for a fixed signal-to-noise ratio, the source SNR roughly evolves as the ratio m/n and therefore improves with m. Second, the manifold-based aggregation may behave quite differently when the dimensionality of the ambient space changes : for some fixed entry-wise error on the estimated mixing matrices, the angular error (i.e. that is proportional to the distance used in the tangent space of the hypersphere) decreases when the number of observations increases.

Figure VI

.1 shows the evolution of the mixing matrix criterion as a function of the mini-batch size for m = 5 (which corresponds to the determined case, as there are n = 5 sources) and m = 20. The sparsity level is fixed to p = 0.1 and the condition number of the mixing matrix A * is equal to 3. The SNR is fixed to 40 dB. At first glance, it is interesting to observe that the dGMCA -Fréchet mean algorithm provides in both experiments decent results when the mini-batch size t b is lower than 80 : here, using small mini-batches enables to alleviate the computational burden while not deteriorating too much the results over GMCA. Being more specific, dGMCA slightly benefits from a larger number of observations m when the mini-batch size is small. This is very likely the consequence of lower errors as measured by the angular distance. Furthermore, for small t b some mini-batches are likely to have a too small number of statistics to be well estimated when m is small. To that regards, the results of the robust Fréchet mean algorithm are particularly informative : while for the smallest mini-batches and small m the results yielded by the mere Fréchet mean were deteriorated, it is not anymore the case. This probably means that the deterioration previously observed with small t b was likely to stem from a small number of badly estimated mini-batches, that is outliers. The use of a robust aggregation lowers their impact, enabling an improvement of about 1 order of magnitude for very small mini-batch sizes (t b < 20). Unexpectedly, equipped with this aggregation procedure, the dGMCA algorithm further allows to improve the separation process with a gain of up to about a factor 2 with respect to the standard GMCA algorithm. We will propose an explanation of such a phenomenon in Section C.3. Lastly, the reader might wonder why dGMCA does not seem to reach the accuracy of GMCA for the largest mini-batch sizes. To explain this, let us recall that the thresholding strategy differs from GMCA to dGMCA : in the case of GMCA, it adapts to the actual distribution of the estimated sources using a percentile. Since such an automatic regularization parameter choice was not applicable in dGMCA, we rather chose a fixed deterministic strategy which is likely to be less effective. When the mini-batch size t b becomes a relatively large fraction of the full number of samples t, the strategy used on GMCA should be also implemented in dGMCA (note that in this setting, it is assumed that the computational burden is not a main issue).

Condition number

The ill-conditioning of the mixing matrix A * , measured by its condition number, plays an important role in the difficulty of a given BSS problem. Mixing matrices with large condition numbers will lead to two major challenges : i) an increased noise level in the source domain 9 , and ii) mixtures that are more colinear. In these experiments, the noise level is fixed to 40 dB and the sparsity level is p = 0.1. A general comment is that as expected the separation quality of all the 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. VI.2). 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 and leads to a gain of about one order of magnitude. This was to be expected, as higher condition number lead to more diversity and less stability among the mini-batches. Similarly, when the mini-batch size decreases, the discrepancy between the two methods increases. Furthermore, while the GMCA results are the best ones when A * is close to orthogonality, the robust dGMCA tends to deteriorate slower as a function of the condition number (cf. Fig VI .3) : that is, here again the aggregation is able to reject the outliers that impede GMCA use (such outliers probably comes from a specific bad noise realization in some mini-batches, which is amplified by high condition numbers). As such, the gain of robust dGMCA over GMCA is very significant for conditions number of about 5. Beyond 10, the curves tends to merge to give bad results, as most of the mini-batches are badly estimated.

Number of sources n

As seen in Chapter V, the number of sources is one of the elements that drive the complexity of blind source separation, as higher n values imply more difficult unmixing problems. In these experiments, the sparsity level is fixed to 0.1, the number of observations For a large number of sources, a general comment is that GMCA obtains degraded performances, consistently with Chapter V. Then, a striking element is that equipped with the Fréchet mean the dGMCA algorithm no more exhibits the peak of improved performances for t b around 20 -50, and eventually does not reach the quality of GMCA in the setting. In contrast, making use of the robust aggregation allows preserving very good separation results for t b < 100, with a gain of almost one order of magnitude with respect to the other methods. These results were to be expected : with a fixed number of observations, finding an increasing number of sources n is more and more difficult, which makes that some mini-batches are strongly badly unmixed. Furthermore, an increasing number of sources also implies an increasing number of partial correlations (i.e. coefficients for which multiple sources with aggregation on Riemannian manifold are active), which are notoriously difficult to handle and make most sparse BSS algorithm deteriorate [Bobin et al. 2015]. From an alternative point of view, a larger number of sources generally increases the number of potential spurious local critical points of the cost function, making robust dGMCA more successful. A more detailed discussion concerning this aspect will be developed in Section C.3.

Sparsity level p

In these experiments, the impact of the sparsity level of the sources p is assessed. Such an impact can be twofold :

-The sources are statistically stationary. However, when they are very sparse (e.g. p = 0.1), the values taken by a single small mini-batch may largely change between two realizations, which can lead to outliers 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 hard-to-escape local minima, which are further smoothed out when sparsity level decreases. ), the GMCA algorithm yields better results and there is far less discrepancy between the two dGMCA me-thods. This highlights that equipped with the (robust) Fréchet mean, combining mini-batch optimization and aggregation leads no performance loss : parallelization allows to go faster and large-scale without diminishing the separation quality. This can be explained by the fact that in this regime the realization of the sources is smoother among mini-batches.

iii) For still less sparse sources (cf. right panel of Fig. VI.5), the results of robust dGMCA are consistent with the previous regime : since the mini-batches are on average still more similar, robust dGMCA obtains a separation quality closer to GMCA for still smaller mini-batches. One of the differences with regime ii) is that GMCA obtains slightly worse results. This is expected, as more partial correlations occur. Furthermore, the discrepancy between Fréchet mean and its robust counterpart increases again strongly. This might also be due to an increased number of partial correlations : a few small mini-batches containing such samples might be much more difficult to unmix.

However, such arguments do not fully answer several observations for the very sparse regime. In particular, why are the results of GMCA enhanced when going from p = 0.1 to p = 0.25 ? Furthermore, why is the robust dGMCA better working than GMCA when p = 0.1 ? Indeed, it should be the opposite, as very large mini-batches should offer in a very sparse regime more stationarity of the precise realization of the sources. To answer to these questions, we will first confirm in the next subsection such observations with a realistic experiment in which the sources are very sparse. In Section C.3, we will then propose an explanation. In brief, the hypothesis we propose is that highly sparse sources induce a less smooth optimization landscape (which is further backed by C.2 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 sensors. In this context, an observation is formed by the linear combination of the contributions from various radionuclides. Each one is described by a signature in energy which is composed of one or several emission lines to which a Compton continuum is associated, as displayed in Figure VI.7. 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 | [END_REF]] and featured in Figure VI.7. The number of observations is fixed to m = 20 and the number of samples per source is equal to 16240. The goal of this section is to evaluate the performances of dGMCA in a setting where the samples are non-stationary, highly sparse and with a large dynamic range ; the information content of a single signature basically spans 2 to 3 orders of magnitude. The sources are modeled in the wavelet domain : γ-ray observations are first decomposed into an undecimated unidimensional wavelet frame [Starck et al. 2010] before applying any BSS method. The number of scales is fixed to 5, which yields a number of wavelet coefficients equal to 81200 ; these are obviously not large-scale data but it already allows to highlight some remarkable results. Figure VI.8 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 solid 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 VI.9 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 very 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 p, 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 hard to escape critical points. 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 along with ALS clearly yields empirical robustness to the algorithm as testified by the very small scatter of the results (the shaded areas). As such, it is likely that the proposed minimization scheme generate some implicit regularization that is beneficial to efficiently tackle sparse BSS problems by making the optimization easier. More concerning this topic is said in the following.

C.3 Discussion -robustness and implicit regularization

During the last few years, understanding the impact of optimization in matrix factorization problems [Gunasekar et al. 2017] or learning (deep) neural networks [Neyshabur 2017] has attracted a lot of interest. More specifically, it has been emphasized in many works that a specific optimization method might enable implicit regularization. In this subsection, we first detail the notion of implicit regularization and then highlight links to such works in order to better understand the behavior of dGMCA. Our hypothesis is backed through further experiments.

Optimization landscape and implicit regularization

Most of the investigations studying the optimization landscape of learning problems focus on the underdetermined or over-parameterized case [Gunasekar et al. 2017, Neyshabur et al. 2017]. In this regime, it has been showed that under some conditions local minimizers are likely to be global. However, such a claim does not help concerning the quality of these critical points [START_REF] Neyshabur | Geometry of Optimization and Implicit Regularizationin Deep Learning[END_REF]. In sparse BSS, the factorization problem at play is furthermore determined or over-determined, which with aggregation on Riemannian manifold larization, either being it explicit or implicit, is of paramount importance. By explicit regularization, we mean a regularization appearing in the cost function (e.g. the sparsity promoting term of Eq. (II.8)). By implicit, we denote an existing regularization only due to the optimization scheme, that is which does not appear explicitly. To give an example of implicit regularization in general matrix factorization problems, one can cite the role induced by gradient descent in many algorithms [Gunasekar et al. 2017]. To tell a long story short, using GD for underdetermined problems implies a regularization penalizing the complexity of the learnt model (e.g. with low nuclear norm)10 . In our case, while such regularization might appear because of the analogy made between ALS and GD in Section B.3, we rather study more specifically the implicit regularization implied by the use of mini-batches.

Regularization induced by mini-batches in machine learning

To that respect, the case of highly sparse sources and more particularly the γspectroscopy example is particularly illustrative : dGMCA performs very well as the batch-size decreases, while GMCA does not. 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 can indeed improve the results over full batch methods [Le-Cun et al. 2012, Keskar et al. 2016, Hardt et al. 2015]. 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 minimizers. Very closely related to dGMCA, several articles have recently emphasized the interest of using non-vanishing learning rates [START_REF] Smith | [END_REF], Xing et al. 2018]. with aggregation on Riemannian manifold This is appealing, as following our analogy with ALS we are specifically in a regime where such a rate is constant over the iterations (set to η (l) = 1). More specifically, and beyond a potential computation gain in terms of the number of parameter updates [START_REF] Smith | [END_REF], it has been empirically shown by [START_REF] Xing | [END_REF]] that using stochastic mini-batches (respectively high learning rate) enables to explore broader areas (respectively to jump optimization landscape borders and escape bad areas), thus finding results further from the initialization. The authors argue that the stochasticity introduced by small size mini-batches through a structured noise strongly favors flat minimizers that are akin to generalize better [Hochreiter & Schmidhuber 1997], since the introduced noise would be roughly aligned with the sharpest loss directions11 .

Links to dGMCA

A similar phenomenon 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, the algorithm performs much better in the highly sparse case, where spurious minima tend to be sharper and more difficult to escape. The stochasticity induced by mini-batches then helps exploring the optimization landscape.

To further investigate such an exploratory ability within the experimental setting of Section C.2, Figure VI.11 displays the histograms of the mixing matrix criterion after 250 iterations of dGMCA (close to "convergence") when the Fréchet mean (left panel) or its robust version (right panel) is used. In the first case, it is interesting to notice that using large mini-batches tends to 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 bad (close to 10 dB), the aggregated estimate is rather poor. Switching to robust aggregation (right panel Figure VI.11), the results might seem to go against the previous hypothesis of exploratory power : while in accordance to the previous experiments the separation with small mini-batches is good, the scatter is however smaller than when using large mini-batches. Nevertheless, we are here looking at the results after 250 iterations. That is, it seems that at this point dGMCA already found a good minimum of the optimization landscape. As such :

-The fact that the scatter with small t b values is smaller that with high t b values indicates that most of the mini-batches are contained within a basin of attraction of a good minimum. This is relieving, as it means that the noise introduced by small mini-batches (and also the high learning rate η (l) ) is not high enough to make most of the mini-batches "jump out" of the good mini-mum. As such, the few outliers present in the histogram are easily handled by the robust aggregation (there is much less than 50 % of them) ;

-The fact that on the other hand the scatter is non-zero, nor even negligible might further highlight that the local minimum we are looking for is rather flat. this is somehow unrealistic in applications since the limiting factor is the total number of samples). The best quartile of the large mini-batches obtains in general almost as good performances as the small ones. Therefore, the lack of exploratory power with large B seems to impede the separation because most of the large mini-batches are unable to move from bad minima. However, some of them are still able to find good solutions, but they are too few for the aggregation to highlight them.

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.

C.4 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.

Complexity of one epoch

Consistently with Chapter V, each iteration of the GMCA algorithm has a complexity of O(t(mn + n 2 + m)). The complexity of one epoch of dGMCA is similar, with aggregation on Riemannian manifold where the last term correspond to the Fréchet mean and L f 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 (see Fig. VI.12) : in practice the computation time for a given number of iterations does not deviate much from linearity (in particular, the transfer costs between the nodes are negligible in comparison to the computation time).

Data reduction ratio Computational time gain

C.5 Experimental results : summing-up

The previous results concerning the computation time have validated the fact that the dGMCA enables to work with large-scale datasets, as there is a linear gain in the complexity of each iteration, and the number of iterations can further be highly reduced in some experimental settings. It was however expected that switching to a distributed version of the GMCA algorithm would not have come at no cost : using mini-batches should have led to larger estimation errors. Seen from the distance, the experiments tend to highlight two distinct regimes : with aggregation on Riemannian manifold -Mildly sparse sources : in this setting, a key result is that 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. This shows that the proposed dGMCA is an useful approach for the separation of large-scale mildly sparse signals, which describe numerous natural signals represented in multiscale representations.

-Sparse to very sparse sources : mini-batch optimization with very sparse sources cause a series of problems. In particular, small-size mini-batches tend to be less homogeneous, which is likely to lead to more outliers and degrade the aggregation procedure. In this case, the robust Fréchet mean is an interesting aggregation procedure. Unexpectedly, it can further continue to improve for very small mini-batch sizes, the best solutions being reached for sizes of the order of the number of sources. Strikingly, the results are then more accurate than the ones of the GMCA algorithm. The explanation we proposed is that dGMCA enables a better exploration of the optimization landscape : if the randomness introduced by the mini-batches is well handled by the aggregation, it makes that the final unmixings are better.

Conclusion

To tackle the large-scale sparse BSS problem, we introduced in this work the dGMCA algorithm, which uses mini-batches in a projected Alternating Least-Square framework. The different estimations of the mixing matrix are aggregated taking into account the Riemannian manifold structure of the problem constraints. To perform such an aggregation, we proposed to use a robust Fréchet mean, which further improve the separation quality. While for mildly sparse sources the approach is experimentally demonstrated to yield a huge gain in computation time at almost no cost in terms of separation accuracy, dGMCA can even improve the separation results over its full batch counterpart when used for very sparse sources. The explanation we propose is that using stochastic mini-batches enables a better exploration of the optimization landscape. Many numerical experiments are proposed to show the relevance of the approach.

Sparse BSS : from Linear to Non-Linear Mixtures

In this last chapter, we extend the previous work by departing from the usual linear setting and tackling the case in which the sources are mixed by an unknown non-linear function. We propose a stacked sparse BSS method enabling a sequential decomposition of the data through a linear-by-part approximation. In this context, non-linear BSS can be seen has solving a potentially large number of linear BSS subproblems. Therefore, an automatic hyper-parameter choice for each sub-problem is mandatory due to the computational burden (cf. Chapter IV1 ), and optimization strategies such as block-coordinate methods (see blockGMCA in Chapter V) or mini-batches (see dGMCA in Chapter VI)) could yield improved performances. Beyond separating non-linearly mixed sources and despite increased indeterminacies compared to linear BSS (cf. Chapter II-A.2.2), the introduced StackedAMCA can under discussed conditions further learn an approximation of the inverse of the unknown non-linear mixing, enabling to reconstruct the sources. The quality of the method is experimentally demonstrated, and a comparison is performed with stateof-the art non-linear BSS algorithms. We also propose an in-depth discussion of StackedAMCA required hypotheses.

A Non-Linear BSS

A.1 Mixing model

In this chapter, the BSS model will change into the non-linear one presented in chapter II-A.2.2 :

X = f * (S * ) + N (VII.1)
Where f * is an unknown non-linear function from R n×t to R m×t (here, n ≤ m). We will here consider general functions f * , by mostly assuming that f * is invertible and symmetrical around the origin, as well as regular enough. Regular means that f * is L-Lipschitz with L small enough and that f * does not deviate from a linear mixing too fast as a function of the input amplitude (which can for instance be the case with sensor saturations, or chemical sensors etc.) ; see Sparse BSS : from Linear to Non-Linear Mixtures on the mixing and the sources, are discussed in Sec. E.

A.2 Previous work

In this subsection, we present some previous works concerning non-linear BSS. A summary of the different algorithm families is displayed in Figure VII.2.

A.3 Independent component analysis

While most of the previous work has been devoted to ICA, the independence prior is not sufficient to ensure source separation in the general non-linear setting [Comon & Jutten 2010]. A first possibility to bypass this separability issue is to explicitly focus on a special kind of mixing f * , for which separability can be shown (under conditions that are not discussed here for the sake of compactness). Among the most well-known kind of mixings, one can cite [Deville & Duarte 2015] :

-Post Non-Linear -PNL. In these, each observed mixture X i , i ∈ [1, m] is an univariate nonlinear function of a linear mixture of the sources :

X i = f * i   n j=1 A * j i S * j   (VII.2)
-Linear Quadractic mixtures -LQ. In these, each observed mixture X i , i ∈ [1, m] is a polynomial function of second degree of the sources :

X i = n j=1 A * j i S * j + n j=1 n k=j B * j i S * j S * k (VII.3)
Where B * is, similarly to A * , a mixing matrix. Another possibility is to use an explicit or implicit regularization making the problem better posed. For explicit regularization, one can cite additional priors on the sources such as temporal dependencies [START_REF] Hyvarinen | [END_REF], Ehsandoust et al. 2017]. However, in this work we will not assume any such additional explicit prior. In such situations, several algorithms use an implicit regularization. Among the most well-known, one can cite [Almeida 2003,Honkela et al. 2007[START_REF] Brakel | Learning Independent Features with Adversarial Nets for Non-linear ICA[END_REF]. The work of [Almeida 2003] is an extension of INFOMAX to the case of non-linear mixtures, in which a neural network is learnt by minimizing the mutual information. The difference with the linear case is that the functions enabling to approximate the cumulative distribution functions of the sources have to be learnt accurately and not chosen a priori, since the non-linear problem is much less constrained. [Honkela et al. 2007] summarizes several works using a variational Bayesian approach bringing the required additional regularization. The derived cost function is used to train a non-linear model of the mixing f * (e.g. a Multi-Layer Perceptron). In [START_REF] Brakel | Learning Independent Features with Adversarial Nets for Non-linear ICA[END_REF], the authors use a Generative Adversarial Network : the idea is to maximize independency measures implicitly thanks to adversarial objectives. The goal of the generator is

B.2 Overview of The Proposed Approach

More specifically, the lowest amplitude samples of the data X can be well approximated by a classical linear model (cf. Fig. VII.3,right) because of the regularity assumption on f * stating that the mixing must not deviate from linearity too fast as a function of the amplitude (cf. Sec. E). Thus, it is possible to find a first linear model for the low-amplitude samples using a sparse linear BSS algorithm, provided that this one is robust to the higher amplitude non-linearities. However, this approximation is too rough for higher amplitude samples, for which the mixing deviates strongly from linearity. It is nevertheless possible to improve it through the introduction of another linear model fitted to higher amplitude samples, thus creating a linear-by-part approximation (with initially two segments). This is while being insensitive to the higher amplitude samples that are more affected by the non-linearities. Indeed, such highly non-linear samples are detrimental to most usual linear sparse BSS algorithms, as i) they behave as partially correlated samples (i.e. samples for which multiple sources are simultaneously active), while many linear BSS algorithms assume disjoint supports of the sources ; ii) they correspond to high amplitude samples, which hinder the use of the Morphological Diversity hypothesis ( [Starck et al. 2010], cf. Chapter III-C.4.2), stating that the highest amplitude samples are the most discriminative for BSS.

R (1) 1 R (1) 1 (a) ( Â(1) † R (1) ) 1 ( Â(1) † R (1)
Due to these two issues, we propose to use the Adaptive Morphological Component Analysis (AMCA - [Bobin et al. 2015], cf. Chap. III-C.5) algorithm, which enables a linear separation of sources having samples with both partial correlations and large amplitudes. In the case of non-linear mixings, the weights M in AMCA enable to discard the samples with high amplitudes which are the most affected by the non-linearities, because these are considered as partial correlations. At iteration l, the algorithm is thus able to find a good linear model Â(l) of the lowest amplitude samples of R (l) . Sparse BSS : from Linear to Non-Linear Mixtures

For selecting only such contributions, there are two issues : i) Determine which samples are well separated by the current linear-by-part model ;

ii) Actually compute a new residual R (l+1) containing the currently badly separated samples, which is done through shrinking their amplitude towards zero.

Solving these two issues is however non-trivial. We could have a fully-sequential naive approach, in which we would re-use at iteration l the previous residual R (l) to compute the new one R (l+1) (this naive approach is illustrated in Fig. VII.5). Issue i) would then be solved by inverting the current linear model found by AMCA, which would align the lowest amplitude samples of R (l) with the axes (cf. Fig. VII.9, where the inversion is performed in the right plot). The badly estimated samples would then be the ones far from the axes, making their selection quite straightforward (cf. The difficulty with such a naive approach however arises when attempting to solve issue ii). Indeed, for each sample the amplitude shrinkage must be performed in the direction of the axis corresponding to the source associated to the sample manifold. For low and middle amplitude samples, this is not a major difficulty : since these are supposed to be well or at least decently unmixed respectively, we can quite reliably assume to which source they belong, and thus in the direction of which axis to shrink. For high amplitude samples, things are however completely different as we have almost no clue concerning the unmixing. Consequently, these can be shrunk in parallel to wrong axes. While this is not a direct issue in the following iteration l +1, during which the linear BSS step will aim only at unmixing the lowest amplitude samples of R (l+1) , this will definitely impact the algorithm through error propagation after several iterations, as R (l+1) will be used to compute R (l+2) , R (l+3) ... In Figure VII.10, we illustrate such a shrinking issue with an example. To determine in which direction each sample must be shrinked, we use an angular criterion (θ = π/4) to associate the samples with the closest axis. In the left plot, the samples in red are associated with the horizontal axis and the samples in yellow with the vertical axis. The shrinking is then operated in the corresponding direction. As can be seen in Figure VII.10 right, the subsequent residual is good for low amplitude samples but bad for higher amplitudes, for which the shrinkage has been performed in the direction of the wrong axis. This will create error propagation in a few iterations, when the algorithm will try to fit a linear model on these high amplitude samples.

B.5.2 Selection Function : proposed solution for computing R (l+1)

To bypass the error propagation issue that would appear in the previous naive fully sequential approach, we introduce at each iteration a preliminary step, that is first presented. The whole solution to issue i) and ii) is then detailed. potential outliers or samples that would be too much affected by noise .

Solution to problem ii) X (l) u can then be shrunk using the values τ (l) 1..n to obtain the residual R (l+1) . For high amplitude samples, as evoked above the current unmixing is bad and we thus do not know in which direction to threshold and which threshold τ

(l) i , i ∈ [1.
.n] to use. Therefore, for these we need to resort to a first guess based on their angular distance to each axis (similarly to Fig. VII.10). While this potentially introduces temporary errors, these are in practice strongly limited due to the two advantages of using the unfolded data X (l) u instead of merely working sequentially and re-using the previous residual R (l) .

B.6 Enhancements

In this subsection, we shortly detail several major enhancements of Stacke-dAMCA that makes it much more robust in practice and enabling to work on more difficult datasets X.

-Correction of the permutations between the linear models Â(1..l) :

The issue is here due to the permutation indeterminacy in usual linear BSS. In our approach, multiple permutations between the linear models Â(1..l) could create a remixing between the sources if they were not properly corrected.

In StackedAMCA, the linear model Â(l-1) found in the previous iteration is used to initialize AMCA at iteration l. Thus, if the mixing f * is regular enough, one can hope that the permutations will not be an issue as the output of AMCA will resemble the input. However, we can introduce a much better permutation correction through re-using the unrolled X (l) u . To do so, we compare after the AMCA step the angle with the axes of the samples that are well linearly separated (i.e. the ones in the red square of Fig. VII.6) to the angle of the corresponding samples with the axes in the unrolled manifold of the previous iteration X (l-1) u

. The permutation correction is then done by minimizing the discrepancy. Note that such a correction is only made possible due to the regularity hypotheses detailed in Section E. It has furthermore been extensively experimentally tested, showing good practical results as confirmed in Section D.

-Determining the sign of the samples on the manifold :

In Sec. B.5, we explained how the direction in which to threshold was determined. A simple approach to determine the way (e.g. from right to left for the samples in the upper right quadrant of Fig. VII.10 but from left to right for the ones in the lower left quadrant) in which to shrink would be to use the sample signs. However, this approach would lead to errors for mixing highly deviating from linearity : in these, the manifolds could cross the axes, making the source sample signs in R n to differ from the ones on the manifold. Therefore, in practice we compute the sign on the manifold of each sample : the process is similar to the one for attributing a source to each sample (cf.

144 Chapitre VII.

Sparse BSS : from Linear to Non-Linear Mixtures solution to problem ii) in the residual computation), but in addition we also attribute a sign.

-Handling non-negative data : Several extensions of NMF to non-linear mixtures have been proposed. In these, the sources are assumed to have non-negative coefficients. Examples of applications are hyperspectral imaging [START_REF] Meganem | [END_REF]

Â(l) , Ŝ 1 2 Tr[(R (l) -Â(l) Ŝ)M(R (l) -Â(l) Ŝ) T ] + n i=1 λ i Ŝi 1 + ι Âi 2 =1 ( Â) +ι . 0 ( Ŝ) (VII.5)
Where the last term enforces an elementwise non-negativity of the coefficients of Ŝ. The rest of the algorithm is the same, except that we only focus on the positive branches of the manifolds when determining the sign of the samples on the manifold4 .

-Enabling different regularities for each manifold :

In the current form of the algorithm, all the manifolds must have the same regularity. That is, their curvature radius must be similar. Otherwise, the most regular manifolds would need less linear models to be well estimated than the other ones. Thus, in the last iterations less sources would need to be separated as the whole most regular manifolds would already have been handled (cf. Fig. VII.12 for a concrete example). We propose to bypass this issue by limiting at each iteration l the amplitude of the thresholds τ (l) 1..n and rather use τ

(l) 1..n = min i∈[1,n] (τ (l) i ) (cf.

C Neural Network Interpretation

StackedAMCA can be interpreted as the multilayer neural network of Fig. VII.13. The different network layers correspond to the iterations l and each layer computes a linear approximation of some part of the data. Therefore, the horizontal manifold is fully handled during the first iteration, while the vertical manifold would need at least two linear models to be well estimated.

(1) u1 X (1) u2 (a) R (2) 1 R (2) 2 (b) X (1) u1 X (
Middle : Residual after the first iteration : only one manifold is left, making the search for two sources irrelevant ; Right : Proposed solution for the threshold choice, to be compared with the left plot.

Slightly altering the notations and writing W (l) = Â(l) † , each neuron layer corresponds to the estimate Â(l) yielded by the linear BSS step. In Fig. VII.13, the non-linear step corresponds to the residual computation. Due to the thresholding, this step has similarities with classical non-linearities in neural networks such as the Rectified Linear Unit (ReLU - [Maas et al. 2013]). This similarity is in particular stronger for the case of sparse non-negative sources described in Sec. B.6, for which a ReLU is truly used to compute the residual. When no non-negativity is enforced, the soft-thresholding operator used resembles a symmetrical ReLU. The network thus roughly possesses the classical alternating between neuron layers and non-linearities. We further need skip connections to complete the transcription of the algorithm. In particular, these enable to re-use X directly, reducing the error propagations and improving the results similarly as usual skip connections [Huang et al. 2017]).

However, contrary to many learning processes the layers are trained one-by-one, each of them minimizing the cost function of AMCA. This would correspond to a greedy training. A global refinement step could however be added. Furthermore, while the non-linear step encompasses a classical non-linearity, it first requires an unfolding of the manifolds, which is different from usual neural networks. Another difference is also that the thresholds τ

D Experiments

D.1 Metrics

Before demonstrating the relevance of the method with three experiments, we present the metrics used to assess the result quality, as these are slightly less usual than in the linear case. Due to the indeterminacy by the non-linear function h in nonlinear BSS, it is important to differentiate measures about the reconstruction of the sources and their separation. The reconstruction metrics are the usual linear ones : the Mean Error (ME), or the more accurate SDR (capturing both the separation and reconstruction quality) [START_REF] Vincent | [END_REF], or SAR (focusing only on the reconstruction). Concerning the separation quality, a classical approach is to estimate h [START_REF] Ehsandoust | [END_REF]] by fitting a non-linear curve P to the 1D-manifold of the scatter plot of each estimated source Ŝi as a function of the true one S * i , and to look at the thickness of the manifold around P. Here, we use for P a polynomial function of degree 20 and the thickness is measured by the :

-logarithmic median absolute distance to P :

C med = -10 log( n i=1 median j (| Ŝj i -P(S * j i ))|) (VII.6)
-logarithmic mean absolute distance to P :

C mean = -10 log( n i=1 1 t t j=1 | Ŝj i -P(S * j i )|) (VII.7)
However, the results of these metrics are sensitive to the choice of P. We thus propose to introduce for exactly sparse sources a new metric based on the angular distance to the axes :

C ang = -10 log    1 n(n -1) n i=1    n i =1 i =i 1 - 1 #Z j∈Z S * j i Ŝj 2 i + S * j 2 i       (VII.8)
where Z = {j ∈ [1, t]|S j i = 0}. Lastly, for the first linear-by-part experiment, it is also possible to use the mixing criterion C A [Bobin et al. 2008], as true underlying linear models exist.

D.2 Linear-By-Part Mixing

The objective of the first experiment is to study StackedAMCA main mechanisms and test it in a ideal setting where it should be able to separate and reconstruct the sources well. The mixing is indeed linear-by-part, thus perfectly matching the unmixing process of StackedAMCA. More specifically, they have t = 10000 samples, with disjoints support and a sparsity level of p = 10%. There is m = n = 2 observations, obtained with a linear-by-part f * (no noise N is added). Each part is an orthogonal A * (l) matrix. The data X is shown in Fig. VII.14 (while the mixing might seem simplistic, it however deviates much from the linearity). Since both the true mixing matrices A * (1..l) and the optimal thresholds are known, it is possible to assess the quality of their estimation by StackedAMCA. Qualitatively and as shown in Figure VII.15 -left, the data reconstruction is almost perfect (except for 2 outliers -likely to come from a thresholding in the wrong direction). This however does not guarantee the separation of the sources. Fig. VII.15-right therefore also displays the scatter plot of X with colors corresponding to the estimated different sources : each manifold is correctly labeled with only one source. Furthermore, the thresholds τ (1..l) 1..n in violet seem to be well estimated (the cumulated error is 4.9 × 10 -3 ). Quantitively, Fig. VII.15 displays the evolution of C A and the SDR as a function of the iterations. The high values of C A confirms that the separation is good, while the decent SDR shows that despite the non-linear setting, the source reconstruction is correct. The high quality results for the first iteration indicates as expected that AMCA is robust enough to discard the highly non-linear high amplitude samples. The decrease of both C A and the SDR with the iteration number is not strictly monotonic, probably due to the fact that some errors done at a given layer l can be compensated at the following layer (e.g. by still finding the good thresholds τ (l+1) 1..n ).

D.3 Star Mixing and Comparison to Other Methods

We now compare the results of our algorithm to other existing ones on a much more complicated experiment. Only a few algorithms for non-linear BSS are open source, and we mostly found three of them : MISEP [Almeida 2003], NFA [Honkela et al. 2007] and ANICA [START_REF] Brakel | Learning Independent Features with Adversarial Nets for Non-linear ICA[END_REF]. The experiment itself is inspired from [Ehsandoust et al. 2016] but is made much more difficult due to a) the presence of noise N with a SNR = 30 dB ; b) a high number of sources (the original experiment being restricted to n = 2). The sources follow a Bernoulli-Uniform distribution, p = 10% of the t = 9500 samples being non-zeros and drawn according to a uniform distribution in [-0.5, 0.5]. We further ensure that the supports of the n = 6 sources are disjoint and there are m = 6 observations. These are computed recursively through the application 15 times on each element j ∈ [1, t] of a mixing of the form U The different separation and reconstruction metrics are displayed in Table VII.1 and VII.2. The corresponding results are shown in Fig. VII.17, where the scatter plot of one estimated source is drawn as a function of the true one for each method. First, it seems that neither ANICA nor NFA truly separate the sources (concerning ANICA, the results seem however to improve when no noise is added). This bad separation is translated into plots that do not resemble a 1D-manifold. It is in particular visible by looking at the estimated values for S * 1 = 0 : many samples are such that Ŝ1 = 0, which corresponds to high interferences and a leakage of the other sources in Ŝ1 . Although such bad unmixings could come from our lack of familiarity with the parameter tuning of these methods, it is possible that the regularization introduced by the network structure for ANICA and the Bayesian setting for NFA is not sufficient to enable the separation of the sources (since the independence is not either, cf. Sec A). On the contrary, MISEP separates the sources well, probably due to a good implicit regularization. StackedAMCA displays a very small number of outliers (a single one in Fig. VII.17). While these probably come from small remaining error propagations due to the temporary first guess in the residual computation (cf. Sec. B), their small numbers shows the interest of the manifold unfolding from the raw data X at each iteration. Therefore, looking at the various metrics the sources are still in general much better separated by StackedAMCA : C mean improves by almost 7 dB compared to MISEP. The good separation of Sta-ckedAMCA is confirmed by the best C ang . Second, MISEP does not reconstruct well the sources as StackedAMCA does and Fig. VII.17 clearly indicates that it did not invert the non-linearity h. On the contrary, the scatter plot yielded by stackedAMCA resembles the identity and the good ME (and best SAR, despite the outliers) indicates that the algorithm structure was sufficient to regularize well the reconstruction problem. Some non-linearities f *

j d = cos(α(j))U j d - sin(α(j))U j f and U j f = sin(α(j))U j d +cos(α(j))U j f , with α(j) = π 2 (1-U j 2 d + U j 2 f ), for random d, f ∈ [1, 6]

E.1 Sparsity of the sources and disjoint supports

The sparsity of the sources S * in the direct domain is assumed for regularizing problem (II.3). Such an assumption can for instance be verified in realistic experiments such as spectrometry.

We have furthermore assumed the supports of the sources to be disjoint. While this is not very realistic in practical cases, it seems difficult to bypass this condition as we only explore the span of f * that the 1D-manifolds created by the sparse sources uncover. By the morphological diversity assumption, the points outside these manifolds are too rare to enable a proper estimation of f * without either further conditions on the mixing (e.g. the separability over the different sources) or additional priors on the sources (cf. Fig. VII.20). Similarly to the results of [START_REF] Hyvarinen | [END_REF], a promising prior could be to use potential temporal dependencies, which is left for future work.

S * 1 S * 2 (a) X 1 X 2 (b) S * 1 S * 2 (c) X 1 X 2 (d)
Figure VII.20 -Illustration of the difficulty to estimate the partially correlated samples in non-linear BSS. Note that with more than n = 2 sources, it is not even possible to determine which sources are active. a) Non-sparse sources : the samples in red are the only ones to correspond to disjoint supports of the sources. All the other ones are partially correlated ; b) A dataset corresponding to the mixing of the sources in a). Since the mixing is not separable, the manifolds in red are not sufficient to estimate the non-linear mixing for the samples in green (the green and red manifolds are not parallel). The non-linearity thus would need to be estimated using the information yielded by the green samples ; c) Sparse sources : the partially correlated samples are in much smaller number ; d) The number of green samples is now too small to estimate the mixing far from the disjoint support samples.

156 Chapitre VII.

Sparse BSS : from Linear to Non-Linear Mixtures

We however did some tests without disjoint supports. The samples with multiple active sources were badly separated but the estimation of the 1D-manifolds by StackedAMCA was not much perturbated, which is due to AMCA robustness to multiple active sources. For instance, re-doing the experiment of Section D.3 without enforcing disjoint supports in the data creation, we obtain similar results : C med = 42 dB, C mean = 26 dB, C ang = 37 dB and -10log(ME) = 27 dB. Note that these metrics were computed on the disjoint samples only since we do not claim to handle well the joint ones, but only that these do not perturbate much the unmixing.

E.2 Density of the 1D-manifolds

For StackedAMCA to work, the manifolds must be continuous, which is to be related with the clustering nature of the algorithm. In particular, this entails that the sources must be dense enough within their support, to avoid any "gap" in the manifolds. The size of affordable gaps depends on the manifold regularity. Furthermore, for each linear BSS step a sufficient number of samples must be non-zero in order to have enough statistics for solving each linear sub-problem (cf. [START_REF] Gribonval | Dictionary Identification-Sparse Matrix-Factorization via 1 -Minimization[END_REF][START_REF] Gribonval | Sparse and spurious : dictionary learning with noise and outliers[END_REF] for more detailed conditions concerning sparse linear BSS identifiability).

E.3 Symmetry of f * around the origin

StackedAMCA currently needs the symmetry of the mixing around the origin. Indeed, after the non-linear shrinkage step, the residual must be almost linear around zero. This might not be verified for non-symmetrical mixings (cf. Fig. VII.21 -the residual is not interpretable as a mixture with n = 2 sources). Such an assumption could in principle be leveraged. First, the data can be symmetrical around a different point as long as a preprocessing step is introduced to center it. Then, tackling nonsymmetrical data could probably be dealt with by introducing non symmetrical non-linear steps (i.e. using a soft-thresholding function with two different thresholds for the positive and negative parts) and using the non-negativity constraint in the linear BSS step. As such, the positive and negative parts of each manifold could be treated separately, and aggregated at the end.

E.4 Well-conditioned linear sub-problems

Each linear sub-problem to be solved by StackedAMCA must be well-conditioned. That is, the corresponding underlying linear model A * (l) must have a high enough condition number so that AMCA can find a good unmixing, which is to be linked with the classical performances of the algorithm in the linear case [Bobin et al. 2015]. Note that however the re-use of the linear model of the previous iteration l -1 as an initialization of AMCA can help it though re-using past knowledge to find a good linear model at the current iteration. An extreme case of badly conditioned (and even colinear) sub-problem at the second E. Discussion : required Hypotheses for StackedAMCA and possible enhancements 157 (2) . The non-symmetrical manifold cannot be tackled as a single source, as it is not linear in zero. The mixing function f * is assumed to be instantaneous, injective, and to depend on the sample amplitude (which is for instance the case with LQ mixtures [Hosseini & Deville 2003, Duarte et al. 2015, Deville & Duarte 2015]). We have furthermore assumed that f * does not deviate too fast from linearity as a function of the amplitude. For differentiable curves, it mathematically means that at every point of the 1D-manifolds described by the mixing X, the local curvature radius must be large enough. This condition is of primary importance to enable StackedAMCA to separate the sources and alleviate the issue of potential permutations between 158 Chapitre VII.

X 1 X 2 (a) R (2) 1 R (
Sparse BSS : from Linear to Non-Linear Mixtures layers (the permutation correction during the iterations is furthermore based on this assumption). For a similar reason, f * must also be L-Lipschitz with L small enough.

E.6 Same length for all the manifolds

All the manifolds must currently have the same length. Otherwise, one could be fully estimated in less iterations than the others. Thus, the residual would comprehend less sources in the last iterations than in the first ones. A possible solution would be to re-estimate at each iteration the number of sources. While the estimation of such a number is a difficult question in BSS, the problem would here be made easier as we would only need to test a smaller number at each iteration. More specifically, at each iteration StackedAMCA could try reduce the number of sources to be estimated by checking if there are enough remaining samples for each manifold in the residual R (l) to perform the separation.

E.7 Low noise

It is important to emphasize that in general, most non-linear BSS algorithm do not assume the mixing to be contaminated by noise (to the notable exception of [Honkela et al. 2007]). While in our approach the noise must be relatively low (cf. [START_REF] Gribonval | Sparse and spurious : dictionary learning with noise and outliers[END_REF] for conditions on the noise concerning the separability of each linear sub-problems), our algorithm seems to be empirically quite robust as demonstrated in Section D. This is mainly due the use of AMCA, which is itself fairly robust [Bobin et al. 2015].

E.8 What Sources can StackedAMCA Reconstruct Well ?

In contrast to other methods, it is further possible to characterize at least one non-trivial type of mixings for which StackedAMCA approximately reconstructs the true sources up to a simple scaling and permutation indeterminacy. A sufficient condition (in addition to the ones required for separability) is that for each sample of the mixing indexed by j, the mixing f * can be written as a product of an unitary matrix and the sources : G = {f * : R n×t → R m×t ∀j ∈ [1, t], X j = A * (S j ) S * j , A * (S j ) ∈ O} (VII.9)

where O is the oblique set. The function A * (.) : R n → R m×n is potentially nonlinearly depending on S * . Since in AMCA the scale of the matrices Â(l) is fixed to 1, if f * ∈ G there is no ambiguity left for the scale of each layer. Due to the regularity assumption, it will then be possible to backproject linearly for each layer the manifold on the axes with small errors and get an approximate reconstruction.

F Conclusion

We have introduced in this Chapter StackedAMCA, a new algorithm tackling the sparse non-linear BSS problem. Based on a new stacked sparse BSS approach, F. Conclusion 159 this method enables to sequentially compute a linear-by-part approximation of the underlying non-linearities. Each linear part is estimated by a robust linear BSS algorithm step, which is followed by a non-linear step . The non-linear step enables to work on increasingly higher non-linearities and is itself composed of an unfolding of the source 1D-manifolds and then a thresholding. We show the relevance of Stacke-dAMCA compared to other state-of-art methods. Beyond separating the sources, in some experiments the algorithm is also able to reconstruct them well despite a severely ill-posed problem. A discussion of the required hypotheses for StackedAMCA to work is furthermore proposed, as well as a characterization of some datasets for which it should be able to reconstruct the sources well.

Chapitre VIII

Conclusion and perspectives

Although sparse Blind Source Separation has well established over the last two decades its capacity to extract meaningful information from multivalued data, most algorithms do not handle properly large-scale problems. The objective of this PhD was to tackle such an issue.

In Chapter IV, we aimed at introducing a robust sparse BSS method using modern-art optimization frameworks such as PALM. While the difficulties to be tackled by such an approach were first discussed, the proposed method enables an automatic choice of the regularization parameters and an increased reliability with regards to the initialization, potentially obviating any relaunch of PALM that would be precluded in the large-scale setting. Chapter V handled mixings comprehending a high number of sources. This problem is especially challenging, as many usual algorithms suffer from decreased separation qualities in such a setting. The introduced bGMCA, based on intermediate-size block coordinate optimization, enables to maintain high quality results while reducing the computational cost. The case of large-scale datasets is dealt with in Chapter VI, in which we introduce mini-batches in GMCA, improving its scalability both in terms of computation time and memory burden. Strikingly, such a method can also in some settings improve the separation quality over full batch ones. Lastly, Chapter VII is an extension of the previous work to the non-linear BSS problem. The proposed StackedAMCA method constructs a linear-by-part approximation of the underlying non-linearities. Numerical experiments highlight its capacity to separate the sources, as well as to reconstruct them in some settings.

The methods have been extensively tested, both on simulated and realistic experiments. As such, they are shown to work well in a wide range of settings.

To foster reproducible research, the codes will be made available online at http: //www.cosmostat.org/software/gmcalab.

Perspectives

The work presented in this thesis have raised several questions, some of which are still open and might be the subject of future research. Conclusion and perspectives

Algorithmic framework and regularization parameter choice

The work presented in Chapter. IV could be prolongated. Although the proposed 2-step approach was shown to work well in many settings, several pathways could be explored :

-Study other optimization schemes. While the use of BCD is shortly discussed in the Appendix E, more work could be dedicated to determine potential differences with PALM, which might appear in more difficult settings than the studied ones (e.g. with Poisson noise, more sources...).

-Use alternative regularization parameter choices than the MAD. In particular, a proper extension of the Stein Unbiased Risk Estimator (SURE) to the sparse BSS problem could be envisioned, as it has led to good results in other inverse problems [Eldar 2009, Giryes et al. 2011]. Preliminary results (in which an already existing SURE method is applied without any further adaptation to the sparse BSS problem) are shown in Appendix E ;

-Beyond the regularization parameter choice, use machine learning tools to enhance the separation.

Use of intermediate block-sizes in sparse BSS

The introduction of bGMCA leads to several questions :

-For the moment, the block choice at each iteration is still relatively unexplored. While we proposed three strategies in Chapter V, they would deserve more work. Finding a good (or even optimal) bloc size is still an open question ;

-Giving more mathematical grounding to the approach. In particular, links with recent works on dropout in matrix factorization [Cavazza et al. 2017, Mianjy et al. 2018] can be highlighted. In these, it is shown that using dropout introduces an implicit regularization promoting low-rank solutions. It has however to be emphasized that bGMCA is quite different from the aforementioned works (in particular, due to the fact that it works on a residual and it is not based on gradient descent but rather on least-squares) ;

-Estimation of the number of sources : such an issue might be difficult in BSS. However, re-using bGMCA and the above remark, it might be possible to perform such an estimation (first results have been obtained in this direction, but are still yet too preliminary to be discussed and must be confirmed) ;

Use of mini-batches in sparse BSS

The use of mini-batches in dGMCA has shown surprising results, making this approach particularly appealing and calling for further extensions :

-Handle datasets with a large number of observations, that is a high number of lines. This is reminiscent in the context of dictionary learning of the recent ODL extension [Mensch et al. 2018]. The main question would then be how to aggregate the different estimations of the sources. The performances of using mini-batches both for the lines and columns would also need to be studied ;

-Combine mini-batches and blocks and study the performances of the resulting algorithm. In particular, this might call for a good block choice, as all the sources might not be present in a given mini-batch ;

-Introduce a refinement stage, akin to what is done the 2-step approach. This should be relatively straightforward as an asynchronous version of PALM has already been proposed in [Davies 2004].

Non-linear sparse BSS

In Chapter VII, we already highlighted some limitations of the current approach and some corresponding possible enhancements. To summarize and add a few elements :

-Building on the geometrical interpretation, we derived required hypotheses for StackedAMCA to work. It would be attractive to mathematically demonstrate these conditions to be sufficient for the separation and the reconstruction (or to find new ones that are) ;

-While similarly to most other works numerical experiments were here conducted on simulations, it would be interesting to apply StackedAMCA on concrete real-life problems ;

-A major issue of the current approach is its inability to handle partially correlated samples properly. To alleviate such an issue, we could use, if they are verified, additional explicit priors on the sources such as temporal dependencies. The median over all the sources Ŝi , i ∈ [1, n] can then be taken to obtain a global criterion. As a side remark, we would like to highlight that such metrics can be highly sensitive to a few badly estimated samples due to the use of a squared 2 -norm.

Non-linear BSS

Metrics for source separation

A classical approach is to estimate the indeterminacy function h [START_REF] Ehsandoust | [END_REF]] by fitting a non-linear curve P to the 1D-manifold of the scatter plot of each estimated source Ŝi as a function of the true one S * i , and to look at the thickness of the manifold around P. The thickness can then be measured by the :

-logarithmic median absolute distance to P : where Z = {j ∈ [1, t]|S j i = 0}.

Metrics for source reconstruction

The reconstruction metrics are the usual linear ones, once the scale indeterminacy corrected : the SDR, SAR... In addition, we also used the Median absolute Error (ME) : Annexe C

ME( Ŝ) =

Convergence conditions

We here detail the conditions ensuring that the algorithms we use converge, if it is the case. Since this work is not new and is not the focus of this thesis, such details were omitted into the main text.

A BCD

Problems of the form of Eq. II.8 can be shown to (almost 1 ) follow the conditions of Lemma 3.1 a) and Theorem 4.1 b) of [Tseng 2001] (cf. example 6.4 therein) :

-The domain of h(A, S) = 1 2 X -AS 2 F is open, and h is differentiable on it. -Writing f (A, S) = 1 2 X -AS 2 F + R S (SΦ T S ) 1 + ι {∀i∈[1,n]; A i 2 2 =1} (A) and taking feasible points 2 for the initialization ( Â(0) , Ŝ(0) ), the level set U = {(A, S); h(A, S) ≤ h( Â(0) , Ŝ(0) )} is compact, and h is continuous on it.

f (A, S) = f (A, S 1 1 , S 1 2 , ..., S t n ) is convex in (S 1 1 , S 1 2 , ..., S t n ) and convex in A. The function f is also regular on U.

-The updates are essentially cyclic, ensuring that each variable is updated enough time. More specifically, in our case the updates are cyclic.

B PALM

It can be proved that our cost function follows the convergence condition of [START_REF] Bolte | [END_REF]].

-The function f defined above follows the Kurdyka-Lojasiewicz property [START_REF] Bolte | [END_REF]. In brief, such property mainly implies that f is sharp (in particular, close to the critical points) up to a reparameterization of its values (there exists a function φ such that the sub-gradients of u → φ • (f (u) -f (ū)) have a norm greater than 1, no matter how u is close to the critical point ū). More specifically, all the terms of Eq. II.8 are semi-algebraic functions -the data 1. To be more precise, the use of ι {∀i∈[1,n]; A i 2 2 =1} (A) makes that the required conditions do not hold, since in particular f is no more convex in A. This is however not an important issue, since the problem only stems for the specific point where A is a matrix having one (or more) columns filled with 0, which is not likely to occur in our context. For more mathematical accuracy, the hyper-sphere constraint ι {∀i∈[1,n]; A i 2 2 =1} (A) can be replaced by the ball constraint ι {∀i∈[1,n]; A i 2 2 ≤1} (A). 2. That is, the initialization respects the conditions enforced by the indicator function.

Annexe C. Convergence conditions fidelity term is a polynomial function, . p is semi-algebraic whenever p is a rational number and the indicator functions are semi-algebraic, making that the sum is semi-algebraic and thus follows the Kurdyka-Lojasiewicz property.

-Both the sparsity and the oblique constraints are proximable.

-The data fidelity term h is C 2 . As such, its gradients ∇h with respect to A and S are Lipschitz continuous on bounded subsets. Furthermore, the corresponding Lipschitz constants are bounded (and non-zeros).

-We use a cyclic update of the variables.

C pALS

The pALS algorithm cannot be proved to converge and can even diverge. In this section, we aim at giving an understanding of such an issue with a simple (historical) example related to NMF. Let us assume that we want to solve the problem, with A ∈ R m×n and x ∈ R m :

argmin s∈R n 1 2 x -As 2 F + ι s∈K + (s) (C.1)
The pALS solution is :

ŝ = Π K + [A † x] = Π K + [(A T A) -1 A T x], (C.2)
which can be, as explained in [Kim et al. 2008], rewritten as :

ŝ = Π K + [ŝ -(A T A) -1 (A T Aŝ -A T x)] = Π K + [ŝ -1 × H -1 (A T Aŝ -A T x)] (C.3)
with H the Hessian of the data fidelity term. Thus, instead of performing a gradient step as it would be the case in PALM, the pALS algorithm performs a quasi-Newton step with projection, but with an arbitrary step 1. Such an update can however increase the cost function as displayed in Fig. C.1, which precludes any convergence guarantee. As a side remark, it is however intuitively visible that the weaker the enforced constraints, the easier for pALS to be able to converge. For instance, when confronted to a sparse minimization problem of the following form :

argmin s∈R n 1 2 x -As 2 F + R S S 1 , (C.4)
pALS is more likely to converge when R S coefficients are small (and in the limit, if R S = 0, it converges). This explains why in practice GMCA tends to stabilize in most cases thanks to its decreasing hyper-parameter strategy, since the hyperparameters are low during the last iterations.

Annexe E

Other exploratory ways of performing a 2 steps strategy

In this appendix, we shortly discuss some alternatives for the refinement step of the 2-step strategy of Chapter IV. This must not be understood as a fully accomplished work, but rather as preliminary results that could pave the way for future research.

A Discussion about BCD

In sparse NMF, the use of a two-step strategy using as refinement a BCD algorithm instead of PALM was proposed in [Rapin 2014]. While using BCD can lead to a computational overload [Chenot 2017[START_REF] Xu | [END_REF] which would not be suitable in our large-scale setting, we tried it on relatively small problems to determine whether this approach could lead to a potential improvement of the separation quality. This section is structured similarly as Chapter IV : we first study the results of an isolated BCD with fixed parameters, and then its behavior as a refinement stage equipped with the MAD heuristic within a two-step approach.

A.1 Fixed parameters

In this section, we use a similar approach as in Chapter IV-C to assess the separation quality of a single BCD. Such a task is again performed using a grid search on the regularization parameters for the data of case 1 and case 2. The protocol is exactly the same as the one we used for PALM. Figure E.1, which is to be compared with Fig. IV.1, indicates very similar results between BCD and PALM and does not seem to suggest that BCD has a better efficiency than PALM when coupled with a grid-search regularization parameter choice : that is, the hyperparameter choice is not less sensitive when using BCD. As a conclusion, such results seem to clearly indicate that a BCD suffers from the same difficult regularization parameter choice than a PALM. This was fully expected on such simple experiments, where the local minima seem to be rare (the standard deviation over the initialization of both PALM and BCD is low). Indeed, since both algorithms are minimizing the same cost function, the only source of discrepancy between both must come from such potential spurious minima. Such preliminary conclusions must however be mitigated, and should be confirmed in more complicated settings with realistic sources, for which such spurious minima could appear, as well as with a higher number of sources n. 1

A.2 BCD as a refinement stage

In this section, we investigate the behavior of BCD within a 2-step strategy as a refinement stage. Similarly to PALM, we use both an initialization by a warm-up GMCA stage and the MAD as an automatic regularization parameter choice within the BCD refinement. The results we obtain in Fig. E.4 for the 2-step BCD do not vary much 2 from the ones of the 2-step PALM of Fig. IV.5. As such, using a BCD does not yield an improvement that would justify the higher computational cost. As a side remark, while for some rare experiments and initializations using the MAD directly within a single PALM could yield decent results as testified by the right plot of Figure IV.5, it does not seem to be the case for BCD. Indeed, in BCD the thresholds computed from a bad initialization are re-used for the whole update of the S matrix, making them more important than in PALM, where they are used only for one proximal gradient step and then updated. This seems to highlight that using a good initialization for a BCD based on a MAD strategy for the regularization parameter choice is even more important than in the case of PALM.

B Decreasing thresholds based on the source distribution

Recent improvements of GMCA enforce more strongly decreasing thresholds by basing their choice on the estimated source distribution (for more details, see previous Appendix). A natural extension of the 2-step strategy is then to replace the MAD by such an enhanced threshold choice based on an increasing percentile of the estimated sources. Preliminary to using it inside of a 2-step approach, a test on a single PALM yields decent results (at least, for a large number of maximal iterations, implying very slowly decreasing regularization parameters -cf. Fig. E.5). The robustness with regards to the initialization is however low compared to the one of the two step approach 3 . Furthermore, a single PALM always gives worse results than the two step algorithm (probably in part due to the fact that the lowering of the high thresholds implies high interferences -this seems to be confirmed by the bad results of error bars). Compared to a continuous decrease, the gain is not as high as expected 4 , which is to be explained by the fact that the reweighting already cancels out most of the interferences appearing due to the gradient step. As such, the strategy by step only brings a limited gain (it however seems to robustify the process with respect to the realization of A * ,S * and N, which was expected since the good thresholds found through GMCA solution are kept during more iterations).

D SURE

The Stein Unbiased Risk Estimator (SURE) is an automatic regularization parameter tuning method based on an unbiased estimate of the MSE of a candidate solution Ŝ [Giryes et al. 2011]. The hyper-parameter is chosen as minimizing such an estimation of the MSE. While originally restricted to additional white Gaussian noise [Stein 1981], it has been generalized to models having the form of an exponential family distribution [Eldar 2009], extending its use to much more various inverse problems. We propose in this subsection to test this approach on the sparse BSS problem. Note that SURE is usually used on non-blind problems, that is problems where  = A * is known in advance, which is not our case here. A full extension of SURE methods to the sparse BSS problem, while interesting, is out of the scope of this work. We rather propose to use the SURE method within the 2-step approach. Following [Giryes et al. 2011], the computation of the threshold 5 could then be based on the currently estimated -at iteration l -Â, which would be supposed to be a decent guess of A * . The SURE estimated parameter should then correspond to an approximation of the best corresponding Ŝ(l+1) in terms of MSE. The hope is that such a Ŝ(l+1) in turns corresponds to a good Â(l+1) and that the refinement step will enter a virtuous circle. Unfortunately this strategy, although itself a cheaper iterative approximation [Giryes et al. 2011] of what a classical SURE method would do, is already too expensive for the large-scale setting since it requires to test many regularization parameters at each iteration. Therefore, we approximated it in our experiments by computing the thresholds only once at the beginning of the refinement stage, based on the first update of GMCA solution 6 . The results of such a method are given in Fig. E.8. While this methods gives relatively decent solutions, the results are deteriorated compared to the κ-MAD two-step approach. This can be understood by the fact that the method does not take explicitly into account the estimation of the unmixing matrix  and thus the threshold choice is only based on noise removal. As such, the 4. Note that the discrepancy increases when the decrease of the thresholds is quicker in PALM -that is, when the number of maximal iterations is smaller.

5. Here, the regularization parameters are the same for all the sources. 6. This strategy thus corresponds to the first update of the thresholds in the iterations of [Giryes et al. 2011]. However, instead of using an initialization based on the golden rule, we here use a pseudo-inverse S(0) = Â † GMCA X instead of S(0) = ÂT GMCA X, consistently to GMCA solution and giving much better results. Annexe F

Elements of Riemannian geometry

Let define M some Riemannian manifold equipped with the metric < , >. One can define for every point x 0 its tangent space T M (x 0 ) as displayed Fig. F.1. The geodesic between two points x 0 and x 1 is defined as its shortest path on the manifold, which then turns to generalize the concept of straight lines on manifolds. In this thesis, the manifold M is assumed to be geodesically complete, which means that there always exists a minimal length geodesic between any two points x 0 and x 1 . The length of the geodesic is then the distance d(x 0 , x 1 ) between these two points. This also entails that one can define for any point x 0 ∈ M the exponential map exp x 0 . This function maps any point of the tangent space at x 0 to M. Under some conditions, the exponential map is bijective and invertible ; its inverse is defined as the logarithm map log x 0 (x 1 ) = exp x 0 (x 1 ) -1 . For some smooth function J : M → R (e.g. the square geodesic distance to only name one), one can define uniquely its gradient at any point x 0 : ∇J (x 0 ) ∈ T M (x 0 ).

  ).mad(u) : Median Absolute Deviation of u. MAD(U) : vector, which elements are the mad of each line of U ;-#U : number of elements in set U ;-exp u and log u : exponential and logarithmic map used for the Fréchetmean ;Specific naming-X ∈ R m×t : dataset ;-A ∈ R m×n : mixing matrix ;-S ∈ R n×t : source matrix ;-N ∈ R m×t : noise matrix ;-P : permutation matrix ;-W = A † ; -U (1) , U (2) , U(3) ... : some generic matrices ; -Θ = {A, S} n : number of sources ;

-

  ∆ : stopping criterion in iterative algorithms ; -C A : mixing matrix criterion (see Appendix B) ; -C med , C mean , C angle : metrics for separation quality. See Appendix B ; s = S * -Ŝ : error on sources ; -E : error on the source estimation introduced by the use of blocks ; -R : residual ; ε : small constant (e.g. 10 -3 ) ; Contents α : parameter of a function (i.e. of a generalized Gaussian distribution, sine or cosine, exponential decay...) ; σ i : standard deviation of the noise in the source indexed by i ; -C d (S) : condition number of S * ; -C d : condition number of the mixing matrix A * ; p : sparsity level (in [0,1]) ; k : number of non-zeros coefficients ; r : block size ; -B : number of mini-batches used ; b ∈ [1, B] : index of a mini-batch ; t b : size of the mini-batch indexed by b ; -J b : indices of the columns in the mini-batch indexed by b ; φ(U (1) , U (2) ) : distance between U (1) and U (2) (e.g. on a geodesic) ; -S m : m-dimensional hypersphere ; ν : parameter of Nesterov smoothing technique ; ρ : step size used for Fréchetmean ; f * : non-linear mixing function ; -X u ∈ R n×t : unfolded manifolds ; τ ∈ R + : a threshold ; h : non-linear indeterminacy function appearing in non-linear BSS (note : in the context of linear BSS, h is a scaling factor) ; -P : polynomial function, supposed to be estimate h ; Non-linear Independent Component Analysis ; -SDR, SAR, SIR, SNR : Signal to Distortion-Artifact-Interference-Noise Ratio ; -ME : Mean absolute Error ; -MSE : Mean Square Error ; -LC / MS : Liquid Chromatography Mass Spectroscopy ; -LC / 1 H NMR : Liquid Chromatography -1 H Nuclear Magnetic Resonance ; -i.i.d. : Independent Identically Distributed -1/2/3/.../nD : 1/2/3/.../D-Dimensional ; -ReLU : Rectified Linear Unit ;

  Figure I.1 -Exemple de problème de BSS en astrophysique, correspondant ici à l'étude d'un rémanent de supernova. Haut : Quelques données / mélanges observés par le satellite Chandra. Chaque image correspond au même rémanent, mais capturé à une longueur d'onde différente. Le mélange peut ici être considéré comme linéaire : chaque image est alors aplanie pour devenir une ligne de X * dans (I.1). Milieu : sources physiques à l'origine des données / des mélanges observés : chaque image (inconnue en pratique et à retrouver par BSS) correspond à une émission spécifique provenant du rémanent. Ici, il s'agit de gauche à droite de : l'émission synchrotron, l'émission thermique et l'émission du fer, à deux décalages vers le rouge différent. Dans (I.1), chaque image est une ligne de S * ; Bas : Chaque courbe est le spectre d'une des 4 émissions précédentes. Ceux-ci (chacun étant une colonne de A * ) ne sont pas disjoints : ainsi, la prise d'images X à différentes longueurs d'ondes ne permet pas d'isoler chaque émission parfaitement, d'où le besoin d'une méthode de BSS.

  Figure II.1 -Example of LC/MS data. Figure taken from [Rapin 2014].

( a )

 a Figure II.2 -Example of BSS problem on the LC/MS data : Left : spectrum of the DL-arginine ; Right : time elution of the DL-arginine. Figures taken from [Rapin 2014].

FigureFigure

  Figure II.3 -Example of BSS with the Chandra satellite : Up : some observed data (each flattened image of 128 x 128 pixels corresponds to one row of X) ; Middle : true physical sources (each flattened image corresponds to one row of S * ), corresponding to several kinds of emissions ; Down : true mixing matrix (each curve corresponds to the spectrum of an emission and is a column of A * )

Figure

  Figure II.5 -Examples of : Left : exactly sparse signal : most of the coefficients are 0 ;Right : approximately sparse signal : while most of the coefficients are non-zeros, many of them are close to 0, making that the signal can be well approximated by its highest amplitude coefficients. Thus, if all the coefficients between the two red lines are set to 0, the remaining samples comprehend most of the signal energy, which makes that the approximately sparse signal is decently approximated by the exactly sparse signal of the left figure.

  Figure II.5b. A good example of sparse signals in the direct domain is the one of the LC/MS data of Fig. II.2. B.1.2 Sparsity in a transformed domain However, most of the signals of interest are not (exactly or approximately) sparse in the direct domain, but rather admit a sparse representation s Φ ∈ R 1×T in a trans-formed domain Φ ∈ R t×T (see Fig. II

  Figure II.6 -Example of use of a wavelet transform, here the starlets, to sparsify an image. Here, the image corresponds to a zoom on a noisy synchrotron emission in a supernovae remnant. Upper left : original image s * ; Upper right : image in the starlet domain S Φ (or more exactly, image corresponding to the first detail scale in this decomposition) : much more pixels are close to 0 ; Down left : histogram of the original image ; Down right : histogram of the image in the starlet domain : the amplitudes are much closer to 0 and the histogram is much less flat than the one in the original domain. In the starlet domain, the data can be considered as approximately sparse.

  Fig. V.1 of Chapter V). Therefore, being able to maintain high separation performances is an open challenge.

  Fig. III.1, where the support of the largest wavelet samples are almost disjoint).

  Just before explaining how GMCA draws on the concept of morphological diversity to disentangle the sources, let us turn towards the geometrical interpretation of sparse BSS by starting back from the example of Fig. III.1. If we draw the scatter plot of the source S * 1 wavelet samples as a function of the ones of S * 2 , we will get a star shape (cf. Fig. III.

Figure

  Figure III.1 -Up : Two natural images which are sparse in the same transformed domain. Down : Wavelet coefficients of the images (to obtain such coefficients, the images have been transformed into grayscale and resized). As hoped due to the morphological diversity hypothesis, since the two images are different the support of their most powerful wavelet samples are different (and here, they are almost disjoint).

Figure

  Figure III.2 -Left, up and middle : Two sources S * 1 and S * 2 . Left, down : Scatter plot of the wavelet coefficients of the two sources. The red square corresponds to the 1 ball of minimum radius enclosing all the source wavelet coefficients. Right, up and middle : Two mixings X 1 and X 2 obtained from S *1 and S * 2 . Right, down : Scatter plot of the wavelet coefficients of the two mixings. The green square corresponds to the 1 ball of minimum radius enclosing all the mixing wavelet coefficients. As expected, its radius is larger than the one of the 1 ball found with the source wavelet coefficients, in red.

  Figure IV.1 -Median of C A (dB) for 5 initializations of PALM algorithm as a 2D function of the 2 thresholds corresponding to the n = 2 sources.

  Fig. IV.1 to highlight the lack of versatility implied by a change of the distribution of S * or of the backprojected noise A * T N.

  Fig. IV.1, each of them corresponding to a new random A * , S * and N. To study the versatility, we look at the diversity among this plots by looking at their dynamic : for each (λ 1 , λ 2 ) value, we plot the mixing criterion of the best estimation minus the mixing criterion of the worst estimation. It gives us an idea of the intra-case variability for given parameters. The results are plotted in Fig. IV.3.

  Figure IV.3 -Dynamic of C A (dB) over different random realizations of A * , S * and N.

  Fig. IV.4. The high values, up to more than 13 dB, are to be compared with the best results of Fig. IV.2b, that is slightly less than 18 dB.

  Figure IV.4 -Dynamic (maximum value minus minimum value) of C A (dB) as a 2D function of the thresholds for 5 different initializations of PALM algorithm.

Figure

  Figure IV.5 -Mixing matrix criterion for PALM and the 2-step strategy (cf. Sec. E) on case 1. The results of GMCA are also plotted as a baseline, with a fixed κ = 3, which is a classical value due to the corresponding hypothesis in terms of Gaussian noise removal. Left : the dashed line is the median of C A over the different A * ,S * and N, and the error bars correspond to the quartiles of the criterion over the initialization ; center : the dashed line corresponds to the values of C A for one specific A * ,S * and N, and the error bars to the quartiles of the criterion over the initialization ; right : the dashed line corresponds to the median of C A over the initializations, and the error bars to the quartiles of the criterion over the realizations of A * ,S * and N.

  The experiments are conducted with 5 realizations of A * , S * and N and 10 different initializations. The results are displayed in Fig IV.5 (red curve). The low efficiency is highlighted by two points : i) the low average values of C A , reaching at most 1 dB in the leftmost plot of Fig. IV.5 ; ii) even in the experiments for which the strategy works, it is very difficult to choose a good κ value, as seen with the upper quartile of the rightmost plot. Furthermore, for one specific experiment (see plot in the center of Fig. IV.5) the standard deviation for the different initializations is very high compared to the average value of C A , which denotes a lack of reliability. Finally, the lack of versatility is shown by the high error bars in the rightmost plot. The goal of the following subsection is to explain these empirical results.

Figure

  Figure IV.6 -Realistic Chandra A * matrix with 12 channels.

Figure

  Figure IV.7 -Up : true sources S * ; Middle : sources estimated by the 2-step algorithm (for a mixing with 30 dB of noise) ; Down : residual between upper and middle images.

Figure

  Figure IV.8 -Left : C A as a function of the SNR (in dB), with a condition number of A * of 10. To plot the figure, the mean over the realizations has been used. The dashed line corresponds to the mean over the different A * and the error bars to the corresponding quartiles. Right : Practical influence of the two terms of Eq. (IV.10) in the experiment.

  Fig. IV.10 for extreme cases close to the corners of the plots of Fig. IV.9, for a condition number of C d = 20. Fig. IV.10a displays a residual resembling closely to a Gaussian distribution, explaining the good 9 results yielded by the MAD within the PALM stage in the corresponding upper left corner of Fig. IV.9. On the contrary, Fig. IV.10d shows a fully non-Gaussian residual since the noise is very low and the A * T A * s term is far from Gaussianity (since α is far from 2), thus explaining the deterioration of the 2-step strategy over GMCA in the lower right corner of Fig. IV

  α

Figure

  Figure IV.9 -Mixing matrix criterion C A (dB) as a function of the SNR and the sparsity level, for 6 values of condition number C d . 2 upper rows : results of GMCA warm-up stage ; 2 bottom rows : improvement over GMCA results yielded by the PALM refinement stage equipped with the MAD strategy.

Figure

  Figure V.1 -Evolution of the mixing matrix criterion C A (cf. Appendix B) of four standard BSS algorithms for an increasing n. For comparison, the results of the proposed bGM CA algorithm are presented, showing that its use allows the good results of GMCA for low n (around 160 dB for n = 3) to persist for n < 50 and to stay much better than GMCA for n > 50. The experiment was conducted using exactly sparse sources S * , with 10% non-zero coefficients, the other coefficients having a Gaussian amplitude. The mixing matrix A * was taken orthogonal. Both A * and S * were generated randomly, the experiments being done 25 times and the median used to draw the figure.

  8) : with more sources in I C , the energy of A I C , ÂI C * , S * I C and ε I C increases, yielding bigger error terms (A I C * -ÂI C )S * I C and -ÂI C s I C . Therefore the errors E become higher, deteriorating the results.

  Fig. V.2b displays the empirical results for several condition numbers C d of the A * matrix. There are n = 50 sources generated in the same way as in the previous experiment : with a Bernoulli-C. Explaining the behavior of bGMCA : numerical experiments on simulated data 95 r (logarithmic scale)

Figure V. 2 -

 2 Figure V.2 -Up : mixing matrix criterion as a function of r for different n. Right : mixing matrix criterion as a function of r for different C d .

Figure

  Figure V.3 -Mixing matrix criterion as a function of r for different sparsity degrees.

r

  Figure V.4 -Right : number of iterations in logarithmic scale as a function of r.

  Figure V.5 -Exemple of A * matrix with 8 columns : the four first columns have spaced means, while the last ones are more correlated

Figure

  Figure V.6 -Up : mixing criterion on realistic sources, with and without a nonnegativity constraint. Down : example of a retrieved source, which is almost perfectly superimposed on the true source, therefore showing the quality of the results.

Figure

  Figure V.7 -Mixing criterion on realistic sources, using a non-negative constraint with r = 10.

  Fréchet mean 1: procedure Fréchet mean on the Oblique ensemble O b (m) 2:

Figure

  Figure VI.1 -Evolution of the mixing matrix criterion with respect to the minibatch size for left : m = 5 and right : m = 20. The number of sources is fixed to n = 5.

Figure

  Figure VI.2 -Evolution of the mixing matrix criterion as a function of the minibatch size for two distinct values of the mixing matrix condition number.

  Figure VI.3 -Evolution of mixing matrix criterion as a function of the condition number for mini-batch sizes t b = 10 and t b = 25.

Figure

  Figure VI.4 -Evolution of the mixing matrix criterion as a function of the minibatch size for 15 sources.

Figure

  Figure VI.5 displays the evolution of C A with respect to the mini-batch size for two values of the sparsity level of the mixing matrix : p = 0.25 in the left panel and p = 0.5 in the right panel. The case p = 0.1 is featured in the left panel of Fig. VI.1. These results seem to reveal three distinct regimes : i) For very sparse sources (left panel of Fig. VI.1), the GMCA algorithm yields decent results, but the robust dGMCA algorithm performs significantly better when small mini-batch sizes are used. In this regime, the realizations of the sources are very different for different mini-batches, explaining the discrepancy between a) robust dGMCA and b) GMCA and Fréchet mean dGMCA ; ii) For mildly sparse sources (left panel of Fig. VI.5), the GMCA algorithm yields better results and there is far less discrepancy between the two dGMCA me-

  Fig. VI.6 : the error bars of GMCA are larger for lower p valuesat least when the results are still acceptable), with potentially more spurious critical points. Using mini-batches enables a deeper exploration of such a landscape, and thus to find a better solution.

FigureFigure

  Figure VI.5 -Evolution of the mixing matrix criterion as a function of the sparsity level with p = 0.25 (left panel) and p = 0.5.

Figure

  Figure VI.7 -Experiment in γ-ray spectroscopy : example of a single observation and the contribution of each of the radionuclide sources.

FigureFigure

  Figure VI.8 -Experiment in γ-ray spectroscopy : 2 2 Na radionuclide estimated 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.

Figure

  Figure VI.10 gives a last interesting insight concerning the optimization landscape of Eq. (II.8). It shows the evolution of the mixing matrix criterion of dGMCA with the robust Fréchet mean as a function of the number of mini-batches B (i.e. this is somehow unrealistic in applications since the limiting factor is the total number of samples). The best quartile of the large mini-batches obtains in general almost as good performances as the small ones. Therefore, the lack of exploratory power with large B seems to impede the separation because most of the large mini-batches are unable to move from bad minima. However, some of them are still able to find good solutions, but they are too few for the aggregation to highlight them.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.

Figure

  Figure VI.10 -Evolution of the mixing matrix criterion as a function of the number of mini-batches for B = 10, B = 25 and B = 100.

C

  Figure VI.11 -Histogram of the mixing matrix criterion across mini-batches with left the Fréchet mean and right the robust Fréchet mean and B = 500 and B = 10.

Figure VI. 12 -

 12 Figure VI.12 -Computational time gain of the dGMCA algorithm with respect to the GMCA algorithm as a function of the reduction gain t/t b . The dGMCA algorithm has been run on a PC cluster equipped with 8 Amd CPUs, each one has 6 cores Istanbul Opteron 8431 at 2, 4Ghz.

  Fig. VII.1 for examples of mixings StackedAMCA can or cannot handle. The other required hypotheses, both 130 Chapitre VII.

Figure

  Figure VII.3 -Left : Scatter plot of n = 2 sources. The red color is associated with samples where only source 1 is active. Yellow is associated with source 2 ; Right : A non-linear mixing of n = 2 sparse sources. The dashed arrows correspond to the mixing directions found by a linear model. The colors, corresponding to each source, are displayed for explaining the distortion introduced by the mixing f * of the source samples but are unknown in a blind setting.

Figure

  Figure VII.9 -Illustration of the AMCA step at the first iteration l = 1. Left : Residual R (1) (note that at the first iteration, R (1) = X) ; Right : Inversion of the linear model found by AMCA : the lowest amplitude samples are well aligned with the axes.

B. 5

 5 Non-linear step : computing R (l+1) B.5.1 Issues for R (l+1) computation The goal of the selection function is to create a new dataset R (l+1) containing only the samples that are not well explained by the current linear-by-part Â(1..l) model, thus paving the way for the next iteration. 140 Chapitre VII.

  Fig. VII.6 left : the badly estimated samples are the ones outside the red square).

  Fig. VII.12 -right for the resulting threshold choice). As such, the condition of equal regularity is transformed into an equal length condition (cf. more about this limitation in Sec. B.6).

X

Figure

  Figure VII.12 -Regularity issue. Left : The two manifolds do not have the same regularity, as the horizontal one is almost linear. The thresholds found by Stacke-dAMCA without any enhancement after the first iteration are displayed in magenta. Therefore, the horizontal manifold is fully handled during the first iteration, while the vertical manifold would need at least two linear models to be well estimated. Middle : Residual after the first iteration : only one manifold is left, making the search for two sources irrelevant ; Right : Proposed solution for the threshold choice, to be compared with the left plot.

  Figure VII.13 -StackedAMCA as neural network

  Figure VII.16 -A 3-dimensional projection of a non-linear mixing with n = m = 6 sparse sources.

Figure

  Figure VII.17 -Scatter plot of one estimated source as a function of the true source : a) StackedAMCA ; b) MISEP ; c) NFA ; d) ANICA.

EFigure

  Figure VII.19 -Scatter plot of one estimated source as a function of the true source. In contrast to Fig. VII.17, the scatter plot of the source obtained by Stacked AMCA does not resemble identity : a) StackedAMCA ; b) MISEP ; c) NFA ; d) ANICA.

Figure

  Figure VII.21 -Example of non-symmetric case. Left : Data X. Right : corresponding residual R(2) . The non-symmetrical manifold cannot be tackled as a single source, as it is not linear in zero.

Figure

  Figure VII.22 -Example of badly conditioned case. Left : Data X ; Right : Residual at iteration l = 2. The two manifolds are not anymore distinguishable.

2 -

 2 Annexe B. Performance metrics for Blind Source Separation causes (e.g. distortions). More specifically :Ŝitarget = Π S * i ( Ŝi ) ŝi interf erences = Π S * ( Ŝi ) -Π S * i ( Ŝi ) ŝi noise = Π S * ,N ( Ŝi ) -Π S * ( Ŝi ) ŝi artif acts = Ŝi -Π S * ,N ( Ŝi )This decomposition is then used to derive the following metrics :-Signal-to-Distortion-Ratio : SDR( Ŝi ) = 10 log 10 Ŝi target 2 2 ŝi interf erences +ŝ i noise +ŝ i artif acts 2 Signal-to-Interferences-Ratio : SIR( Ŝi ) = 10 log 10 to-Artifacts-Ratio : SAR( Ŝi ) = 10 log 10 Ŝi target +ŝ i interf erences +ŝ i noise

  absolute distance to P :C mean = -10 log( Ŝj i -P(S * j i )|) (B.4)However, the results of these metrics are sensitive to the choice of P. We thus used for exactly sparse sources a new metric based on the angular distance to the axes :

Figure

  Figure C.1 -A graphical example of an iteration of pALS, during which the cost function does not decrease. The gray ellipses represent the level sets of the datafidelity term. The algorithm starts at the current iterate ŝ(k) . The update before thresholding s reaches the minimum of the data fidelity term. However, the projection of s on the non-negativity constraint increases the whole cost function (which is confounded with the data-fidelity term for feasible points) since the new estimation s(k+1) goes from an inner ellipse towards an outer one. Figure inspired from [Kim et al. 2008].

  Fig. E.2 and Fig. E.3, which are compared with Fig. IV.3 in the right plots, tend to show that the versatility of BCD is further similar as the one of PALM, at least close to the diagonal where the best results are achieved.

Figure

  Figure E.1 -Median of C A (dB) for 5 initializations of BCD algorithm as a 2D function of the 2 thresholds corresponding to the n = 2 sources. To be compared with Fig. IV.1. Left : Case 1, Right : Case 2.

Figure

  Figure E.2 -Dynamic of C A for BCD over different random realizations of A * , S * and N in case 1. To be compared with Fig. IV.3 left plot, in the case of PALM. Left : BCD results ; Right : Absolute difference between left plot and the corresponding one for PALM of Fig. IV.3.

FigureFigure

  Figure E.3 -Dynamic of C A for BCD over different random realizations of A * , S * and N in case 2. To be compared with Fig. IV.3 right plot, in the case of PALM. Left : BCD results ; Right : Absolute difference between left plot and the corresponding one for PALM of Fig. IV.3.

Figure

  Figure E.7 -Results on case 2 of PALM with a percentile-by-step threshold choice. For comparison, the results of PALM with a continuous percentile decrease of the thresholds are displayed, as well as the usual 2-step κ-MAD algorithm. The upper plots correspond to isolated PALMs and the lower ones to PALMs included in a two-step algorithm. Left : the dashed line is the median of C A over the different A * ,S * and N, and the error bars corresponds to the quartiles of the criterion over the initialization ; Right : the dashed line corresponds to the median of C A over the initializations, and the error bars to the quartiles of the criterion over the realizations of A * ,S * and N.

Figure F. 1 -

 1 Figure F.1 -Sketch of some Riemannian manifold M and its tangent space at the point x 0 .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 B PALM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 C pALS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 D Variants of GMCA thresholding strategy 175 A κ-MAD with varying κ . . . . . . . . . . . . . . . . . . . . . . . . . . 175

	vi		Contents
	VIII Conclusion and perspectives	161
	A Proximal operators	165
	B Performance metrics for Blind Source Separation	167
	C Convergence conditions	171
	A	BCD . .

B Increasing percentile . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 E Other exploratory ways of performing a 2 steps strategy 177 A Discussion about BCD . . . . . . . . . . . . . . . . . . . . . . . . . . 177 A.1 Fixed parameters . . . . . . . . . . . . . . . . . . . . . . . . . 177 A.2 BCD as a refinement stage . . . . . . . . . . . . . . . . . . . . 178 B Decreasing thresholds based on the source distribution . . . . . . . . 180 C Decreasing by steps threshold . . . . . . . . . . . . . . . . . . . . . . 181 D SURE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 F Elements of Riemannian geometry 187 Bibliography Contents 1

  Tableau IV.2 -C A for 5 SNR values and 5 algorithms.

		10 dB 15 dB 20dB 30 dB 60 dB
	2 step	15.0	16.3	17.4 19.7	20.9
	PALM 11.9	13.3	13.5	14.2	14.5
	GMCA 13.2	14.8	15.1	17.1	18.6
	EFICA 8.8	10.3	14.0	18.9	19.4
	RNA	9.8	12.6	15.6	18.3	18.4

  ,Eches & Guillaume 2013, Meganem et al. 2014] or show-through removal[Liu & Wang 2013]. A possible extension of StackedAMCA is therefore to handle properly the case of nonlinear sparse non-negative sources (for realistic examples of such sources, see e.g.[Rapin 2014]). Due to the linear-by-part structure of the algorithm, such an enhancement is quite straightforward as it mostly amounts to transform the sparse linear BSS step into a NMF sparse linear one, making it possible to re-use all the corresponding literature (see e.g.[Rapin et al. 2013, Rapin et al. 2014]). Therefore, instead of looking for a minimizer of Eq. (VII.4), AMCA will here aim at minimizing :

	min

  MethodC med C mean C ang slightly worst results compared to the previous subsection, as the problem is more complicated due to the different manifold regularities. StackedAMCA however still obtains the best ones, which are good as testified by Figure VII.19. Therefore, as expected the algorithm is still able to perform the separation.Concerning the reconstruction, it can already be seen in FigureVII.19 that the nonlinearity is not well inverted. This is confirmed by a decrease of 8.4 dB of the SAR and a loss of 3.1 dB for the ME (to be compared with a loss of only 0.8 dB for C mean ).In conclusion, even in this experimental setting where the source reconstruction cannot be expected from StackedAMCA and is indeed bad in practice, the algorithm is still able to separate the sources well.

	StackedAMCA 37.4	27.3	38.1
	MISEP	23.0	22.2	21.2
	NFA	18.4	11.4	12.2
	ANICA	16.7	12.0	4.14

  difficulty of the problem at hand and the fact that most of the previous works based on ICA and general mixings did not provide required conditions under which the source separation can be hoped[Almeida 2003[START_REF] Brakel | Learning Independent Features with Adversarial Nets for Non-linear ICA[END_REF]. In contrast, the geometric interpretation of StackedAMCA enables to perform such a discussion. Concerning the other works based on sparsity, the works making the hypotheses explicit also used strong ones (cf. e.g.[Puigt et al. 2012]). Furthermore, most of the discussed conditions are either intrinsic to the problem or can be mitigated. Hypotheses E.1 and E.2 are intrinsic to sparse non-linear BSS based on clustering (at least, without any additional explicit priors), as well as E.5 (with the regularity potentially depending on the specific clustering algorithm used). While conditions E.3, E.4 and E.6 are specifically required by our algorithm, we propose exploratory paths for alleviating E.3 and E.6. Assumptions E.7 and E.8 must be seen as advantages of StackedAMCA over most other existing algorithms. Concerning E.7, most algorithms do not incorporate any noise in the mixing process. Thus, being able to deal with noise is a progress, even if this one should be limited to some extent. Concerning E.8, to the best of our knowledge StackedAMCA is the first non-linear BSS algorithm for which the source reconstruction of the sources can be hoped for a characterized class of non-trivial mixings.

	154 Chapitre VII.	Sparse BSS : from Linear to Non-Linear Mixtures
	Tableau VII.4 -Reconstruction quality of StackedAMCA, MISEP, NFA and
	ANICA.
	Method	-10log(ME) SAR
	StackedAMCA	23.3	11.0
	MISEP	20.9	7.46
	NFA	13.7	5.51
	ANICA	2.09	-16.9

Où A * (taille m × n) est une matrice contenant les coefficients linéaires de mélange. De manière très simplifiée, l'objectif de la séparation de sources est de démélanger les sources S * à l'origine des données (compte tenu de certaines indéterminations qui seront détaillées dans le reste du manuscrit). Ceci est d'autant plus délicat que la séparation aveugle ne suppose (quasiment) aucune connaissance a priori sur le processus de mélange f * , si ce n'est dans le cas où le mélange est supposé linéaire, et où celui-ci revient donc à une multiplication matricielle. En l'état, le problème de BSS formulé comme précédemment admet une infinité de solutions, parmi lesquelles seul un petit nombre correspond aux signaux physiques, c'est à dire qui sont réellement à l'origine des données observées. Dit autrement, la BSS est un problème mal posé, pour lequel il est courant de rajouter une information a priori dans le but de restreindre l'espace des solutions possibles.1. Les abréviations utilisées dans cette introduction correspondent, par soucis de cohérence avec le reste du manuscrit, aux acronymes anglais.

Exacte est à comprendre ici dans le sens où ces algorithmes garantissent de converger vers un point critique du problème (I.2). En revanche, ce point n'est pas nécessairement un minimum global du problème.

Par heuristique, nous entendons ici une méthode approximative simplifiant un des aspects du problème initial (réduction du temps de calcul, choix d'hyper-paramètres ou d'initialisation...), et permettant d'obtenir des résultats corrects mais non nécessairement optimaux.

-L'introduction d'une méthode permettant de traiter des jeux de données de grandes tailles, et en particulier présentant un grand nombre d'échantillons t. Le problème a ici principalement trait à aux difficultés calculatoires, car un grand nombre d'échantillons i) provoque un temps de calcul fortement accru, rendant très couteuses des méthodes telles GMCA[Bobin et al. 2007] puisque celles-ci nécessitent, à cause de leurs structures itératives, de multiples inversions de matrices gigantesques ; ii) peut impliquer des tailles de données ne permettant même pas de stocker les matrices en jeu en mémoire, prohibant l'utilisation des méthodes de BSS

; 4 -Une extension à la séparation aveugle de sources non-linéaire. Peu de travaux ont été proposés en BSS non-linéaire parcimonieuse. La méthode introduite durant ce doctorat nécessite de résoudre plusieurs sous-problèmes de BSS linéaire, ce qui peut dans des cas complexes potentiellement nécessiter de résoudre auparavant les points précédents.

https ://www.skatelescope.org/the-ska-project/

https ://www.computerworld.com.au/article/392735/ska_telescope_generate_more _data_than_entire_internet_2020/

Note that we only detail the two proximal splitting methods that will be used in this work. Nevertheless, a vast literature related to this topic exists, with many well-known algorithms : Alternating Direction Method of Multipliers (ADMM -[START_REF] Boyd | [END_REF]), Douglas-Rachford[Douglas & Rachford 1956],[START_REF] Chambolle | [END_REF], Dykstra-like splitting[START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]...

Among these extensions, some accelerated versions of PALM have been proposed (cf. e.g.[Pock & Sabach 2016,Hien et al. 2019]). In this work and due to the difficulty of performing sparse BSS using PALM (cf. Chapter IV), we do not however consider such accelerations, which is left for future works.

We point out that while notations such as A(∞) and S(∞) are slightly abusive in this context because GMCA is not proved to converge, the goal is to clarify the rational of hyper-parameter choice.

Due for example to an implicit regularization introduced by the minimization scheme.

While the case of n = 2 sources might seem too simplistic, the following experiments will already highlight the difficulty of performing sparse BSS using PALM in this setting. As such, it is not expected that the use of PALM would be made easier with more sources.

http ://chandra.harvard.edu/

Introducing a weighted Frobenius norm -squared Mahalanobis distance -would avoid this situation but it requires a good estimate of the covariance matrix of the noise N.

Note however that while such a conclusion seems to hold true in this case, where the condition number of A * is relatively low and where the sources are not too sparse, it might not be the case for all experiments, see the plots of next subsection.

National Institute of Advanced Industrial Science and Technology (AIST), Spectral database for organic compounds : http://sdbs.db.aist.go.jp

Depending on the instrumentation, high SNR values can be reached in such an experiment.

When using ODL, the regularization parameter choice has to be performed by the user. In the simulations we propose, we tried several values and kept only the best results.

Note that beyond matrix factorization, similar implicit regularizations have been studied in machine learning and it has been emphasized that in contrast to standard belief GD algorithm leads to minima[Lee et al. 2016] that surprisingly generalize well[Gunasekar et al. 2017,Neyshabur 2017, Gidel et al. 2019, Azizan et al. 2019]

Note that for similar reasons, some authors have been injecting noise in their algorithm to benefit from such a better exploration of the optimization landscape[Neelakantan et al. 2015]see[Rapin 2014] in sparse BSS.

Note that we however only use only the warm-up stage ; the introduction of a refinement step is left for future work.

Note that in this case, the symmetry of the mixing is not anymore required.

While such decent results could suggest to use a two-step PALM -PALM with the first stage using a threshold choice based on the decreasing strategy, and the second one using a MAD strategy in addition to reweighted 1, the results would likely be unstable with regards to the initialization similarly to the first step.

Remerciements

Chapitre IV. Using PALM in Sparse BSS

X i 0 = cos(α(i)) 5 S * i 0 -sin(α(i)) 5 S * i 1 and X i 1 = sin(α(i)) 5 S * i 0 + cos(α(i)) 5 S * i 1 ; Middle : Mixing deviating from linearity, but regular enough for StackedAMCA to work well : X i 0 = cos(α(i))S * i 0 -sin(α(i))S * i 1 and X i 1 = sin(α(i))S * i 0 + cos(α(i))S * i 1 , with α(i) = π 2 (1 -S * i 2 0 + S * i 2 1 ) ; Right : Linear mixing : X = A * S * .

132 Chapitre VII.

Sparse BSS : from Linear to Non-Linear Mixtures Figure VII.2 -Summary of the different families of algorithms. In red, settings we study in this chapter. The "YES" mention for separability might be subject to other assumptions. The "CLAIMED" refers to the work of [Ehsandoust et al. 2016].

then to produce samples reconstructing the data well (therefore, a reconstruction is performed through a decoder) and which are undistinguishable by the discriminator of independent samples. The implicit additional prior thus comes here from the network structure itself.

A.4 Sparsity

Concerning works based on the sparsity of the sources, [START_REF] Theis | [END_REF], Van Vaerenbergh & Santamaría 2006, Duarte et al. 2015] attempted the problem of PNL mixtures, and [Duarte et al. 2012] the problem of LQ mixtures. General settings similar to the framework of the current article have mainly been studied in [Ehsandoust et al. 2016,Puigt et al. 2012] and focus on the geometric interpretation. In [Ehsandoust et al. 2016] the authors claimed the separability of the sources when using a sparsity prior, if only one source is active for each sample. Note that the mere sparsity prior is however in general not sufficient for source reconstruction. The algorithm they propose uses a clustering approach, followed by a manifold unfolding (which, in particular aims at tackling the partially correlated samples if additional information is known -e.g. source regularity). However, they do not propose a specific method to perform the clustering in the case of more than n > 2 sources. In [Puigt et al. 2012], the authors first propose to find zones in the observation domain in which only a single source is active. This gives a scattered representation of the mixing non-linearities. As several zones may correspond to the same source, these are then aggregated to get more statistics for estimating the mixing. Unfortunately, the experiments also only focus on the case n = 2.

A.5 Contribution

We propose to tackle the general problem of non-linear BSS presented in Eq. (II.3) by using a sparsity prior on the sources, without assuming any additional explicit priors on them. To the best of our knowledge, our method is the first attempting to find a linear-by-part approximation of the underlying non-linearities using a stacked sparse BSS approach. This departs from usual methods as :

-In contrast to neural network approaches [Almeida 2003[START_REF] Brakel | Learning Independent Features with Adversarial Nets for Non-linear ICA[END_REF],

we explicitly use the geometric interpretation in terms of 1-dimensional (1D) manifolds (cf. Sec. B) existing in the case of sparse sources ;

-In contrast to clustering approaches [Ehsandoust et al. 2016,Puigt et al. 2012], the clustering we use is only local and linear-by-part. As such, it is easier in settings with a high number of sources : in particular, we will apply it in the case of more than n > 2 sources contrary to [Ehsandoust et al. 2016, Puigt et al. 2012]. Furthermore, it is based on usual sparse linear BSS algorithms, enabling to re-use much of the corresponding work in terms of optimization and automatic hyper-parameter choice, which has lead to enhanced separation quality [Bobin et al. 2015] (cf. previous Chapters).

Beyond separating sources, the algorithm yields a possible source reconstruction by inverting the estimated linear-by-part model. Despite the usual non-linear BSS indeterminacies, this reconstruction is empirically shown to estimate well the true sources under some discussed hypotheses. In Sec. B, the method is further described. In Sec. D, some experiments are conducted on three different mixings, two of them with a high number of sources. In Sec. E, the required hypotheses for the proposed approach are studied.

B Proposed Approach

B.1 A Geometrical Perspective on Sparse Non-Linear BSS

The proposed method is first described by adopting a geometric point of view in the case of n = 2 sources. The principle can however be generalized to higher values, and the relevance of the method will be empirically testified on such difficult settings in Section D. When plotting S * 1 as a function of S * 2 , most of the source samples lie on the axes due to the morphological diversity hypothesis ( [Bobin et al. 2007] -cf. Chap. III-C.4.2). In this chapter, we will even assume that all the source samples lie on the axes (cf. Fig VII.3,left) -this disjoint support hypothesis is further discussed in Sec. E. Once mixed through the non-linear f * function, the source samples on the axes are transformed into n non-linear 1D manifolds [Ehsandoust et al. 2016, Puigt et al. 2012], each corresponding to one source (see Fig VII.3,right).

B.3 StackedAMCA, detailed description and notations

The whole algorithm iterates the process described in the previous subsection by alternating at each iteration l a linear BSS step and a non-linear step. The first step computes a linear model Â(l) on the current residual R (l) . The second one paves the way for the next iterations by computing a new residual R (l+1) . This is done by finding for each source i a maximum amplitude value τ (l) i above which the non-linearities are too high to be considered as currently well estimated, and then shrinking the current data using the τ (l) 1..n (cf. Fig. VII.6 for the first residual computation). This shrinkage enables to sequentially reduce the amplitudes of the originally higher non-linearities, and therefore to compute linear models describing them. Thus, at each iteration a new linear model is fitted to increasingly higher amplitude samples of the original data X.

More details concerning the two main steps are given below. The algorithm as well as the notations are summarized in

B.4 Linear Sparse BSS Step : AMCA

The main requirement for the linear sparse BSS algorithm is its abillity to find a linear model representing well the lowest amplitude samples of the residual R (l) , 136 Chapitre VII.

Sparse BSS : from Linear to Non-Linear Mixtures -Linear step : estimates Â(l) with AMCA through the minimization of (cf.

Sec. III-C.5 for more details) :

u , select highly non-linear samples : find τ showing the interest of such an unfolding to enhance the first guesses about the unmixing. The sample in red can be seen to lie much closer to the good axis after the l = 3 iterations, justifying the simple use of an angular distance to the different axes to make a decent first guess in the residual R (l+1) computation.

Preliminary step : unfolding the manifolds. Instead of working directly on the previous residual, at each iteration l the algorithm starts back from the raw data X and unfolds the manifolds using all the previously computed linear models Â(1..l) and τ (1..l-1) (cf. u . This approach has two advantages 2 : -Less error propagation occurs in the iterative process, as at each iteration we start back from the raw data X, instead of R (l) ; -The contrast between the sources for high amplitude samples is increased, which makes that the separation is easier. For instance, in Fig. VII.11, the red dot is naturally made closer to the good axis after l = 3 iterations through the re-use of all the previously computed linear models.

Solution to problem i) StackedAMCA uses for each sample of X (l)

u the distance (or more specifically the angle) to the axes. If such a distance is small enough, it is assumed that the sample is well separated by the current linear-by-part model. For each source i, we denote τ (l) i the maximum amplitude of the samples close 3 enough to the axis of i. In practice, the choice of the τ (l) i is made much more robust by using a clustering method based on the amplitude of the samples and enabling to discard 2. In addition to these two advantages, the unrolled X (l) u is also used to correct possible permutations between the different linear models Â(1..l) , cf. Section B.6

3. Such a closeness criterion depends on the manifold regularity, and is one of the few algorithm hyper-parameters to be tuned by the user.

D.4 Experiment without source reconstruction

The objective is here to try stackedAMCA on a mixing for which it should not be able to reconstruct the sources well. The experimental setting is relatively similar to the one of the previous subsection : the sources follow a Bernoulli-Uniform distribution, p = 10% of the t = 9 500 samples being non-zeros and drawn according to a uniform distribution in [-0.65, 0.65]. We further ensure that the supports of the n = 6 sources are disjoint and there is m = 6 observations. The mixing is computed recursively through the application 9 times on each element j ∈

, where all the variables where defined in the previous subsection. The most important difference is due to the presence of the exp(α(j))/2 factor, making that the first column of the corresponding mixing matrix is not anymore unitary, and thus stackedAMCA cannot be expected to perform source reconstruction (cf. Sec. E). Furthermore, the presence of this term on the first column only makes that the different manifolds do not have anymore the same regularity, enabling also to test StackedAMCA in this setting. A 3D-projection of the mixing is displayed in Fig.

Proximal operators

Definition of proximal operators

The proximal operator of an extended-valued proper and lower semi-continuous convex function f : R t → (-∞, ∞] is defined as :

Where u ∈ R t .

Definition of the soft thresholding operator

For two scalar λ and b, the soft thresholding operator S λ (.) is defined as :

For two matrices Λ and U, we extend our notation to :

Definition of the projection of the columns of a matrix on the 2 unit hypersphere

We define the projection of a column vector u on the 2 unit hypersphere as :

We extend the notation to a matrix U by projecting all its columns on the 2 unit hypersphere :

Definition of the projection of a matrix on the positive orthant

Annexe B

Performance metrics for Blind Source Separation

The choice of a performance metric is key to quantitively assess the quality of a given separation. In this Appendix, we summarize the criteria used in this thesis. First, we highlight that finding such a criterion might not be trivial, and has lead to several propositions [Bobin et al. 2015, Vincent et al. 2006, Ehsandoust et al. 2017]. For instance, a mere recovery of the data (e.g. X ÂŜ ) does not mean that A * and S * are well estimated. Second, we will categorize the metrics into two families : the ones used for linear BSS, and the one used in the non-linear setting, for which the separation and reconstruction quality are distinct.

Linear BSS

The most largely used metric in this work is the mixing matrix criterion C A [Bobin et al. 2015] :

With A * the true mixing matrix and P Â † the pseudo-inverse of an estimate corrected through P for the scale and permutation indeterminacies. The mean is the average of all the elements inside the matrix. The smaller C A , the better the separation. However, in most of this work we used for the sake of clarity a logarithmic criterion : -10 log 10 (C A ).

Note that due to the use of the pseudo-inverse of Â, this metric might be sensitive. Furthermore, C A does not allow to distinguish between a fair but not fully accurate estimate of A * and a good estimate in which for instance only one column is badly estimated. As such, an angular criterion between the columns of  and A * might be sometimes preferred [Chenot 2017].

Originally introduced in the context of audio BSS, the SDR, SAR, SIR, SNR [START_REF] Vincent | [END_REF]] are also of interest. In brief, the authors decompose each of the estimated source Ŝi , i ∈ [1, n] in :

Where the four terms should be seen as respectively the part of Ŝi coming from the wanted source Ŝ * i , the one coming from other sources, from noise and from other Annexe D

Variants of GMCA thresholding strategy

We here present a few variants of the originally-proposed GMCA automatic parameter choice based on the use of the κ-MAD. The main idea is the same, and all these strategies aim at accentuating the decrease of the parameters to still further benefit from the morphological diversity assumption.

A κ-MAD with varying κ

In this strategy, the hyper-parameters are chosen as described in chapter III on the κ-MAD rule :

However and contrary to what was described, in this strategy κ is not chosen fixed with a value of κ = 3 but decreases linearly during the iterations, starting for instance from a value of 7 and finishing to a value of 3. This enables to better take into account the morphological diversity, since the hyper-parameters now start from higher values.

B Increasing percentile

The value of κ = 7 of the previous strategy might seem somewhat arbitrary. Therefore, this method is fully based on the estimated S(l) i coefficients : -First the support of S(l) i is determined, keeping only the samples of S(l) i larger than κ MAD( S(l) i ) (with κ = 3) ; -Denoting L the maximum number of iterations, the threshold value µ

Compared to the previous one, this strategy has two advantages : i) it avoids the need for setting κ to an arbitrary value ; ii) it enables to take into account the distribution of the sources. However, it has to be highlighted that since the thresholds values are dependent on the final number of iterations L, the results are in general much more dependent to L than with the previous strategy. a decreasing threshold strategy within PALM with exactly sparse S * sources, for which the decrease is faster). Such a strategy is however made more robust to the initialization when used within a 2-step approach. Indeed, the choice of the initial regularization parameters is then made on a decently estimated Ŝ matrix, making it much better. Since GMCA is robust to the initialization, so do the threshold choice. The final results are however slightly worse (2.5 dB on average) than the ones of the 2-step κ-MAD approach presented in chapter IV. An interesting result is furthermore given when using a 2-step strategy with a refinement stage based on BCD instead of PALM and a threshold choice based on the percentile instead of the MAD. In this setting, BCD yields better results than a refinement stage based on the PALM. This can be explained by the fact that in BCD S is updated until convergence with fixed thresholds before updating A. Thus, the interferences implied by the gradient step are partially removed before updating A, implying better results. While we do not aim at performing a full study of BCD, this result suggests to use within PALM a decreasing by steps threshold choice, aiming at partially mimicking BCD.

C Decreasing by steps threshold

This new strategy is based on the previous remark that a decreasing threshold choice based on an increasing percentile of the estimated sources worked better within a 2-step BCD than a 2-step PALM, and that it might be linked to lower interferences in BCD. More specifically, high thresholds creates high artifacts due to the use of the 1 norm. If the thresholds decrease quickly, such artifacts will be directly partially transformed into interferences due to the PALM gradient step on S. In turns, these interferences will impact the estimation of A, deteriorating it, which is limited by the use of BCD. While this seems to be relatively intrinsic to the PALM scheme, we can at least try to limit the deterioration by improving the threshold choice and not choosing it on a solution containing a high level of interferences. More specifically, instead of decreasing the threshold at each iteration, we will compute a new threshold only after several iterations, when the interferences become low. That is, we will decrease the thresholds by step, by choosing when to change the thresholds using a convergence criterion on Ŝ. The results of this method are displayed in Fig. E.7. Such a strategy yields deteriorated results with a single PALM compared to the previous decreasing strategy. This is understandable, since at the beginning of the algorithm the bad initialization implies bad thresholds, which are kept during more iterations in a decreasing by steps scheme. On the other hand, this strategy gives similar results as the usual κ-MAD rule (while they are in fact slightly better, the difference is not significative with respect to the Abstract : During the last decades, Blind Source Separation (BSS) has become a key analysis tool to study multi-valued data. The objective of this thesis is however to focus on large-scale settings, for which most classical algorithms fail. More specifically, it is subdivided into four sub-problems taking their roots around the large-scale sparse BSS issue: i) introduce a mathematically sound robust sparse BSS algorithm which does not require any relaunch (despite a difficult hyper-parameter choice); ii) introduce a method being able to maintain high quality separations even when a large-number of sources needs to be estimated; iii) make a classical sparse BSS algorithm scalable to large-scale datasets; and iv) an extension to the non-linear sparse BSS problem. The methods we propose are extensively tested on both simulated and realistic experiments to demonstrate their quality. In-depth interpretations of the results are proposed.

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France