Diagonalisation conjointe nonnégative par congruence pour l'analyse en composantes indépendantes semi-nonnégative La Diagonalisation Conjointe par Congruence (DCC) d'un ensemble de matrices apparaît dans nombres de problèmes de Séparation Aveugle de Source (SAS), tels qu'en Analyse en Composantes Indépendantes (ACI). Les développements récents en ACI sous contrainte de nonnégativité de la matrice de mélange, nommée ACI semi-nonnégative, permettent de tirer profit d'une modélisation physique réaliste des phénomènes observés tels qu'en audio, en traitement d'image ou en ingénierie biomédicale. Par conséquent, durant cette thèse, l'objectif principal était non seulement de concevoir et de développer des algorithmes d'ACI semi-nonnégative basés sur de nouvelles méthodes de DCC nonnégative où la matrice de passage recherchée est nonnégative, mais également d'illustrer leur intérêt dans le cadre d'applications pratiques de séparation de sources.

II

Contents 4 Blind source separation applications 4.1 Separation of simulated MRS signal . . . . . . . .

Chapitre 1: Introduction

La SAS consiste à estimer P sources inconnues à partir de N observations. Ces observations sont des mélanges linéaires instantanés bruités des dites sources satisfaisant le modèle suivant [Comon and Jutten, 2010]:

x = As + ε (1) où x = [x 1 , x 2 , . . . ,
x N ] T ∈ R N est le vecteur d'observations, où s = [s 1 , s 2 , . . . , s P ] T ∈ R P est le vecteur de sources, où A ∈ R N ×P est la matrice de mélange et où ε = [ε 1 , ε 2 , . . . , ε N ] T ∈ R N est un vecteur de bruit additif. Le terme "aveuglement" fait référence au fait que les sources et la matrice de mélange sont inconnus. La SAS montre son intérêt dans un grand nombre d'applications, incluant les télécommunications, l'imagerie biomédicale et l'acoustique. Le problème de SAS est mal posé sans hypothèses supplémentaires. Afin d'obtenir des résultats pertinents dans une application réelle, il est nécessaire d'utiliser certains a priori sur les sources ou le modèle de mélange. Exploiter l'indépendance statistique des sources à parti d'un nombre fini de réalisations du vecteur aléatoire x (1) caractérise l'ACI [Jutten andHérault, 1991,Comon, 1994]. L'hypothèse de nonnégativité des sources ou de la matrice de mélange, voire les deux, a donné naissance à un large éventail d'approches de SAS, telles que la Factorisation Matricielle Nonnégative (FMN) [Lee and Seung, 1999] et l'ACI nonnégative [Plumbley, 2003]. Dans ce chapitre, nous avons effectué un bref survol bibliographique des méthodes de SAS les plus connues, à savoir l'ACI, l'ACI nonnégative et la FMN. En outre, la méthode de DCC et son lien avec l'ACI sont été décrits. Ensuite, nous avons introduit le modèle proposé d'ACI semi-nonnégative [Wang et al., 2013,Coloigner et al., 2014a,Coloigner et al., 2014b]. Plus précisément, dans le modèle d'ACI semi-nonnégative, la matrice de mélange (contrairement aux sources) est supposée avoir des composantes nonnégatives. En fait, selon ce modèle, les tranches matricielles des tableaux cumulant suivent le modèle de DCC:

C (k) = AD (k) A T , k ∈ {1, . . . , K} (2)
où A ∈ R N ×P + est la matrice de passage nonnégative et égale à la matrice de mélange, où D (k) ∈ R P ×P sont des matrices diagonales, et où K est le nombre total de tranches matricielles du tableau cumulant. Par conséquent, on peut estimer la matrice A à l'aide de méthode de DCC sous contrainte de nonnégativité sur la matrice de passage, conduisant à la notion de DCC nonnégative. Durant cette thèse, nous avons proposé de nouveaux algorithmes de DCC nonnégative pour résoudre le problème d'ACI semi-nonnégative.

Chapitre 2: Algorithmes de diagonalisation conjointe nonnégative par congruence outre, la complexité numérique pour chaque algorithme a été analysée, ce qui permet une comparaison équitable de leurs performances.

Chapitre 3: Simulations numériques

Les performances des algorithmes proposés de DCC nonnégative ont été évaluées sur des jeux de matrices synthétiques, conjointement diagonalisables. Cinq méthodes de DCC non-orthogonale classiques sans contrainte de nonnégativité, incluant ACDC [Yeredor, 2002], CLU [Maurandi et al., 2013], FFDIAG [Ziehe et al., 2004], LUJ1D [Afsari, 2006], LUJ2D [Afsari, 2006] et une méthode de DCC nonnégative nommée LM + sym [Coloigner et al., 2014c], ont été testées comme méthodes de référence. Les performances ont été évaluées en termes de précision de la matrice estimée et de complexité numérique. Tout au long du chapitre, i) les propriétés de convergence, ii) l'influence des erreurs de modèle, iii) l'incidence du nombre de matrices à diagonaliser, iv) l'effet de la cohérence des vecteurs colonnes de la matrice de passage, et v) l'influence du module d'unicité, ont été largement étudiés au moyen de réalisations de Monte Carlo. Les résultats ont montré qu'une meilleure précision d'estimation pouvait être obtenue en particulier dans des contexte difficiles, en exploitant l'information préalable de nonnégativité. Généralement nos algorithmes de type Jacobi fournissent de meilleurs résultats, pour un problème difficile de DCC nonnégative, par exemple pour de faibles valeurs de rapport signal-sur-bruit, pour un petit nombre de matrices à diagonaliser, pour des niveaux élevés de cohérence de la matrice de passage, et pour une grande valeur du module d'unicité. Parmi eux, l'algorithme basé sur le critère des moindres carrés indirects sans étape de prétraitement offre le meilleur compromis entre précision d'estimation et complexité numérique. Lorsque le problème de DCC nonnégative considéré est bien conditionné, c'est-à-dire lorsqu'une valeur élevée du rapport signal-sur-bruit est considérée ou un grand nombre de matrices à diagonaliser est disponible, l'algorithme basé sur l'optimisation des directions alternées est le plus efficace sur le plan de la précision d'estimation pour une complexité numérique équivalente. En la comparant à une méthode existante de DCC nonnégative basée sur l'optimisation de Levenberg Marquardt [Coloigner et al., 2014c], dont la performance est tout aussi compétitive pour certaines expériences, les méthodes proposées requièrent moins de temps de calcul (environ 3 fois moins). De manière générale, les algorithmes proposés offrent le meilleur compromis entre la performance et la complexité.

Chapitre 4: Applications de séparation aveugle de source Dans ce chapitre, nous avons appliqué nos algorithmes de DCC nonnégative à un ensemble de tranches matricielles d'un tableau cumulant d'ordre quatre, donnant naissance à une nouvelle classe de méthodes d'ACI semi-nonnégative. L'ACI semi-nonnégative se retrouve dans un certain nombre de problèmes de SAS impliquant une matrice de mélange nonnégative, par exemple dans i) l'analyse de composés chimiques en spectroscopie par résonance magnétique, ii) l'identification des profils spectraux d'harmoniques (par exemple, de notes de piano) d'un morceau de musique mono-canal par décomposition du spectrogramme, iii) l'élimination partielle du texte se trouvant au verso d'une feuille de papier fin. Les méthodes proposées ont été appliquées pour résoudre ces problèmes afin de vérifier leur capacité de séparation de sources. Nos méthodes ont été comparées à cinq algorithmes de référence en SAS, à savoir deux méthodes d'ACI: CoM2 [Comon, 1994] et SOBI [Belouchrani et al., 1997], trois méthodes exploitant la contrainte de nonnégativité: ACI nonnégative [Plumbley, 2003], FMN [Kim and Park, 2008] et semi-FMN [Ding et al., 2010]. Au vu des résultats obtenus, les méthodes classiques, surtout l'ACI nonnégative, SOBI et semi-FMN, n'étaient parfois pas capables de séparer les sources latentes. CoM2 et FMN ont pour leur part donné de bons résultats d'estimation, malgré quelques imperfections. Les méthodes proposées se sont comportées de la même manière et ont été plus aptes à fournir de meilleurs résultats que les méthodes classiques, en ce qui concerne la précision d'estimation des sources. Nos méthodes ont une complexité numérique modérée lorsque le nombre d'observations n'est pas trop grand. Cependant, leur complexité croît très vite lorsque le nombre d'observations augmente. Heureusement, en incorporant une étape de prétraitement de compression nonnégative, la complexité globale des méthodes proposées d'ACI semi-nonnégative reste acceptable. Ces applications de SAS ont démontré la validité et l'intérêt des algorithmes proposés. Les résultats expérimentaux ont aussi montré qu'en exploitant pleinement l'information préalable des données, telles que la nonnégativité de la matrice de mélange et l'indépendance statistique des sources, les méthodes proposées fournissent de meilleurs résultats d'estimation.

Chapitre 5: Conclusion et perspectives

Ce dernier chapitre conclu le rapport de thèse offrant une discussion sur les différentes contributions de la thèse, mais aussi sur les possibles extensions à venir.

Nonnegative joint diagonalization by congruence for semi-nonnegative independent component analysis

I

Rèsumé

La Diagonalisation Conjointe par Congruence (DCC) d'un ensemble de matrices apparaît dans nombres de problèmes de traitement du signal, tels qu'en Analyse en Composantes Indépendantes (ACI). Les développements récents en ACI sous contrainte de nonnégativité de la matrice de mélange, nommée ACI semi-nonnégative, permettent de tirer profit d'une modélisation physique réaliste des phénomènes observés tels qu'en audio, en traitement d'image ou en ingénierie biomédicale. Par conséquent, durant cette thèse, l'objectif principal était non seulement de concevoir et développer des algorithmes d'ACI semi-nonnégative basés sur de nouvelles méthodes de DCC nonnégative où la matrice de passage recherchée est nonnégative, mais également d'illustrer leur intérêt dans le cadre d'applications pratiques de séparation de sources. Les algorithmes de DCC nonnégative proposés exploitent respectivement deux stratégies fondamentales d'optimisation. La première famille d'algorithmes comprend cinq méthodes semi-algébriques, reposant sur la méthode de Jacobi. Cette famille prend en compte la nonnégativité par un changement de variable carré, permettant ainsi de se ramener à un problème d'optimisation sans contrainte. L'idée générale de la méthode de Jacobi est de i) factoriser la matrice recherchée comme un produit de matrices élémentaires, chacune n'étant définie que par un seul paramètre, puis ii) d'estimer ces matrices élémentaires l'une après l'autre dans un ordre spécifique. La deuxième famille compte un seul algorithme, qui utilise la méthode des directions alternées. Un tel algorithme est obtenu en minimisant successivement le Lagrangien augmenté par rapport aux variables et aux multiplicateurs. Les résultats expérimentaux sur les matrices simulées montrent un gain en performance des algorithmes proposés par comparaison aux méthodes DCC classiques, qui n'exploitent pas la contrainte de nonnégativité. Il semble que nos méthodes peuvent fournir une meilleure précision d'estimation en particulier dans des contextes difficiles, par exemple, pour de faibles valeurs de rapport signal sur bruit, pour un petit nombre de matrices à diagonaliser et pour des niveaux élevés de cohérence de la matrice de passage. Nous avons ensuite montré l'intérêt de nos approches pour la résolution de problèmes pratiques de séparation aveugle de sources. Pour n'en citer que quelques-uns, nous sommes intéressés à i) l'analyse de composés chimiques en spectroscopie par résonance magnétique, ii) l'identification des profils spectraux des harmoniques (par exemple, de notes de piano) d'un morceau de musique mono-canal par décomposition du spectrogramme, iii) l'élimination partielle du texte se trouvant au verso d'une feuille de papier fin. Ces applications démontrent la validité et l'intérêt de nos algorithmes en comparaison des méthodes classique de séparation aveugle de source.

V

Acronyms and Notations

Acronyms

Signal processing area Notations

Typography • R N 1 ×N 2 ו••×N i denotes the set of real-valued (N 1 × N 2 × • • • × N i ) arrays • R N 1 ×N 2 ו••×N i + denotes the set of nonnegative real-valued (N 1 ×N 2 ו • •×N i ) arrays
• vectors are denoted by bold lowercase letters (a, b, . . .)

• matrices are denoted by bold uppercase letters (A, B, . . .)

• HO arrays are denoted by bold calligraphic letters (A, B, . . .)

• the (i, j)-th entry of a matrix A is symbolized by A i,j

• the (i, j, k)-th component of a third order array B is symbolized by B i,j,k

• A :,i (also a i for simplicity) denotes the i-th column vector of matrix A; A i,: denotes the i-th row vector of matrix • stands for the matrix element-wise division

• T is the transpose operator

• -1 is the inverse operator for square matrices

• -T denotes the inverse after the transpose operator

• is the Moore-Penrose pseudo-inverse operator

• vec(A) reshapes a matrix into a column vector by stacking its columns vertically

• det(A) is the determinate of a square matrix A

• diag(A) returns a matrix comprising only the diagonal elements of matrix A

• Diag(b) gives a diagonal matrix whose diagonal elements are given by the vector b

• off(A) vanishes the diagonal elements of the input matrix A

• E{.} denotes the mathematical expectation

List of Tables

1.1 Some different assumptions of NMF, semi-NMF, ICA, NICA and SeNICA for solving the BSS problem x = As + ε. . . . . . . . . . . . . . . . . . . 1.2 The differences between the proposed six nonnegative JDC methods according to the cost functions, the optimization schemes and the assumption on A. When we use the square change variable A = B 2 , the LU and QR matrix factorizations are used to decompose B, leading to the first five Jacobi-like optimization methods. The sixth method based on ADMM optimization estimates A without change of variable.

[1] Regarding iJDC + LU-1 and iJDC + QR-1 , the joint transformation matrix A is constrained to be square. And a preprocessing step that inverses each input matrix C (k) is required. . . . . . . . . . . . . . . . . . . . . . . . . . . . . [Yeredor, 2002], CLU [Maurandi et al., 2013], FFDIAG [Ziehe et al., 2004], LUJ1D, LUJ2D [Afsari, 2006] and LM + sym [Coloigner et al., 2014c The Blind Source Separation (BSS) problem consists of retrieving unobserved sources from the observations. In this thesis report we focus on the following linear instantaneous mixture model in the field R of real numbers:

Numerical complexities of the proposed methods and of ACDC

x = As + ε (1.1)
where x = [x 1 , x 2 , . . . , x N ] T ∈ R N is the observation vector, where s = [s 1 , s 2 , . . . , s P ] T ∈ R P is the source vector, where A ∈ R N ×P is a mixing matrix, and where ε = [ε 1 , ε 2 , . . ., ε N ] T ∈ R N is an additive noise vector. The blindness refers to the fact that neither the sources nor the mixing matrix is known.

The BSS problem appears in a wide range of application areas, where in general a multi-sensor system is involved, including remote sensing [START_REF] Ma | A signal processing perspective on hyperspectral unmixing: Insights from remote sensing[END_REF], telecommunication [START_REF] Chevalier | Application to telecommunications[END_REF], multiple-input multiple-output radar systems [START_REF] Li | Adaptive blind source separation and equalization for multiple-input/multiple-output systems[END_REF]], medical imaging [Wang et al., 2010], brain-computer interface [START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF], Albera et al., 2008], image analysis [START_REF] Be | Blind separation of superimposed shifted images using parameterized joint diagonalization[END_REF], and audio separation [START_REF] Fuentes | Harmonic adaptive latent component analysis of audio and application to music transcription[END_REF]. For example, in the context of biomedical signal processing, the cutaneous ElectroCardioGraphy (ECG) measured on the mother's skin is a non-invasive technique that enables possibility of visualizing the electrical activity of a fetal heart, namely the fetal ECG. These cutaneous ECG recordings can be considered as an approximate noisy instantaneous linear mixture of potential signals generated by underlying bioelectric phenomena, such as the maternal and fetal heart activity, potential distributions generated by respiration and stomach activity [START_REF] Zarzoso | Fetal ECG extraction from maternal skin electrodes using blind source separation and adaptive noise cancellation techniques[END_REF]. Therefore, in order to evaluate the health condition of the fetus and to reveal information for early diagnosis, BSS methods were proposed to estimate the fetal ECG from the recordings on the mother's skin [START_REF] Zarzoso | Fetal ECG extraction from maternal skin electrodes using blind source separation and adaptive noise cancellation techniques[END_REF],De Lathauwer et al., 2000, Niknazar, 2013].

In the BSS model (1.1), it is easy to see that the scale and permutation indeterminacies will necessarily hold. More precisely, a scaling factor can always be exchanged between one source s i and the corresponding column vector of the mixing matrix A, namely a i , without changing the product. Similarly, permuting the order of the sources and the corresponding column vectors of the mixing matrix will not change the observations. Hence the amplitude, the sign and the order of the sources are not uniquely defined.

Therefore, BSS methods estimate the source vector up to a permutation and a scale diagonal matrix. Besides the scale and permutation indeterminacies, the BSS problem (1.1) is still ill-posed without additional assumptions on the model. In order to achieve relevant separating results with an actual application, it is necessary to utilize prior knowledge about the sources' characteristics or about the mixing model. For example, the sources could satisfy some basic assumptions, such as statistical independence, nonnegativity or sparsity. Similarly, the mixing matrix could also preserve some properties, such as nonnegativity or sparsity. The fundamental statistical independence assumption of the sources leads to the well-established theory of Independent Component Analysis (ICA) [Jutten andHérault, 1991, Comon, 1994]. The nonnegativity assumption on the sources or on the mixing matrix or on both, gives birth to a wide range of successful BSS approaches, such as Nonnegative Matrix Factorization (NMF) [Lee and Seung, 1999] and Nonnegative ICA (NICA) [Plumbley, 2003]. During this thesis, the main objective was to develop BSS methods that retrieve statistical independence sources assuming a nonnegative mixing matrix. That is to say the mixing matrix in model (1.1) is constrained to have nonnegative components, while the values of the sources are unconstrained, giving rise to the concept of Semi-Nonnegative ICA (SeNICA). By now only a few methods have been proposed to address this problem [START_REF] Coloigner | Semi-nonnegative independent component analysis: the (3,4)-SENICA exp method[END_REF], Coloigner, 2012, Coloigner et al., 2014a, Coloigner et al., 2014b]. In this report, we propose to solve the SeNICA problem by means of a class of new Nonnegative Joint Diagonalization by Congruence (NJDC) methods.

The rest of this section is organized into two parts. In the first part we provide a brief overview of some well-known BSS methods. Section 1.2 summarizes the existing results of ICA. Section 1.3 is devoted to get a closer look on Joint Diagonalization by Congruence (JDC) and its connection to ICA. Afterwards we introduce some nonnegativity constrained BSS approaches, including NMF in section 1.4 and NICA in section 1.5. In the second part the main principle of the proposed methods is presented. We address the proposed SeNICA model as well as the problem of NJDC in section 1.6. Furthermore, a list of my publications is given in section 1.7.

Independent component analysis

Generally, the methods based on the assumption that the sources are non-Gaussian (or at most one Gaussian source), and statistically independent, belong to the category of ICA approaches. ICA has emerged into a rapidly growing field due to the pioneering work of Jutten and Hérault [START_REF] Hérault | Space or time adaptive signal processing by neural network models[END_REF], Jutten, 1987, Jutten and Hérault, 1991] and later to the work of Comon [Comon, 1994] who gave a mathematical framework of ICA. The ICA problem is defined as follows:

Problem 1.1. Given M realizations of an observation random vector x ∈ R N , find a mixing matrix A ∈ R N ×P and the M corresponding realizations of a source random vector s ∈ R P , such that:

x = As + ε (1.2)
where s has statistically independent components and ε ∈ R N is an additive noise random vector independent of s.

Chapter 1

Independent component analysis

Statistical independence is a reasonable assumption in many realistic problems, and ICA finds its applications in various areas including, but not limited to, radiocommunications [START_REF] Chevalier | Application to telecommunications[END_REF], biomedical engineering [START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF], Albera et al., 2010,Albera et al., 2012,Hajipour Sardouie et al., 2014], remote sensing [START_REF] Moussaoui | On the decomposition of Mars hyperspectral data by ICA and Bayesian positive source separation[END_REF], and face recognition [START_REF] Bartlett | Face recognition by independent component analysis[END_REF], Kim et al., 2005]. Furthermore, many convenient measures of independence were proposed, such as kurtosis, negentropy, mutual information and cumulants, leading to different ICA algorithms.

In the rest of the section, after a brief discussion about the identifiability and preprocessing step of ICA, we review the following four popular ICA methodologies: ICA by maximization of non-Gaussianity, ICA by maximum likelihood estimation, ICA by minimization of mutual information, and ICA by tensor decomposition of cumulant arrays. For a broader survey we refer to [Comon and Jutten, 2010]. Overviews of some recent advances of ICA can be found in [Hyvärinen, 2013, Adalı et al., 2014].

Identifiability and preprocessing

In fact, if the ICA model (1.2) is noiseless, it is essentially identifiable under the following three conditions [Comon, 1994, Eriksson and[START_REF] Eriksson | [END_REF]:

1. The joint probability density function (pdf) of s factorizes as a product of the marginal pdfs of the sources s i .

2. The pdf's of the source components s i are non-Gaussian.

3. The mixing matrix A must be of full column rank.

The identifiability signifies that the estimated mixing matrix and sources are equal to the original mixing matrix and sources, respectively, up to the permutation and scale indeterminacies. The second assumption could be possibly relaxed for at most one Gaussian component. The identifiability issues of noisy ICA were addressed in [Davies, 2004].

Generally it is assumed that the observation vector has zero mean, which can be done by centering it. Moreover, optionally the use of Principal Component Analysis (PCA) as a spatial whitening has been adopted as a mean to reduce the search space of the mixing matrix to the group of unitary matrices, and also to reduce the effect of additive noise. PCA transforms the observation vector x ∈ R N into a vector z ∈ R P by multiplying it on the left a whitening matrix W ∈ R P ×N , such that W (x -ε) has an identity covariance matrix. The resulting mixing matrix W A ∈ R P ×P is an orthogonal matrix.

The row vectors of W are the P eigenvectors of the covariance matrix of x associated with the P largest magnitude eigenvalues, scaled by the inverse of the square root of the corresponding eigenvalues. It is noteworthy that the whitening step introduces additional errors that cannot be fully compensated for in the following ICA step.

Non-Gaussianity, maximum likelihood and mutual information

In this section, we review three fundamental ICA ideologies, including the ICA algorithms based on non-Gaussianity maximization, the ICA methods based on maximum likelihood estimation, and the ICA approaches based on mutual information minimization.

Chapter 1

ICA by maximization of non-Gaussianity

The central limit theorem says that the distribution of a sum of statistically independent independent and identically distributed (iid) random variables tends toward a Gaussian distribution. The fundamental idea is that, if none of the sources has a Gaussian pdf, the estimated source s j = V j,: x usually has a distribution that is closer to Gaussian than any of the original sources s i and that becomes the least Gaussian when s j is equal to one of the original source s i , i.e. when V j,: is the j-th row vector of the demixing matrix V . Therefore maximizing the non-Gaussianity of s j = V j,: x with respect to V j,: results in finding one of the sources. In practice, the kurtosis and the negentropy are two commonly used measures of non-Gaussianity.

• Kurtosis. The kurtosis of a random variable s, denoted by kurt(s), is defined as follows:

kurt(s) = E{s 4 } -3(E{s 2 }) 2 (1.3)
assuming s having zero mean, where E{•} denotes the mathematical expectation.

The absolute value of kurtosis is zero for a Gaussian variable, and greater than zero for most non-Gaussian random variables [START_REF] Hyvärinen | Independent component analysis: algorithms and applications[END_REF]. Many ICA algorithms were proposed by maximizing the kurtosis. Delfosse and Loubaton proposed a deflation approach [START_REF] Delfosse | Adaptive blind separation of independent sources: A deflation approach[END_REF]. Hyvärinen and Oja gave a fast fixed-point algorithm in [START_REF] Hyvärinen | A fast fixed-point algorithm for independent component analysis[END_REF]]. Papadias proposed a stochastic gradient method based on a multiuser kurtosis criterion [Papadias, 2000]. Zarzoso and Comon introduced a line search optimization for the kurtosis contrast function [START_REF] Zarzoso | Robust independent component analysis by iterative maximization of the kurtosis contrast with algebraic optimal step size[END_REF]].

• Negentropy. The differential entropy function H(•) of a random vector s with a joint pdf denoted by p s is defined as follows [START_REF] Cover | Elements of Information Theory[END_REF]]:

H(s) = -p s (η) log(p s (η))dη (1.4)
Based on equation (1.4), the negentropy function J(•) is defined as follows:

J(s) = H(s Gaussian ) -H(s) (1.5)
where s Gaussian is a Gaussian random vector of the same covariance matrix as s.

Negentropy is always nonnegative and it is zero if and only if s has a Gaussian distribution. The computation of negentropy is difficult since it needs to estimate the pdf. Therefore, approximations of negentropy were proposed in order to simplify the computation, for example [Hyvärinen, 1997]:

J(s) ∝ (E{G(s)} -E{G(ν)}) 2 (1.6)
where G(•) is a non-quadratic functions, and where ν is a vector of zero-mean unit-variance Gaussian variables. Based on the above approximation, Hyvärinen proposed a fast fixed-point ICA algorithm, known as FastICA [Hyvärinen, 1999]. Koldovský et al. provided an improved version of FastICA whose accuracy attained the Cramár-Rao lower bound [START_REF] Koldovský | Efficient variant of algorithm fastICA for independent component analysis attaining the Cramár-Rao lower bound[END_REF].

Chapter 1

1.2. Independent component analysis

ICA by maximum likelihood estimation

For the noiseless mixture x = As, if the pdf is specified for each source, the distribution for the whole data set x is determined for any value of A. It is denoted by p x|A and the Maximum Likelihood (ML) estimate of A is A ML = argmax A p x|A [Cardoso, 2010].

The ML is computationally complicated, since it requires to estimate the pdf's of the sources. There are two ways of avoiding the estimation of pdf's. The first way is to use the prior densities in the likelihood. The second way to is to approximate the densities of the variables using a family of simple densities that are specified by a limited number of parameters as described in [Girolami, 1998, Lee et al., 1999, Hyvärinen et al., 2001].

The ICA approach by ML estimation was originated by Pham et al. [START_REF] Pham | Separation of a mixture of independent sources through a maximum likelihood approach[END_REF], and gained popularity after the introduction of the Information Maximization (InfoMax) principle by [START_REF] Bell | An informationmaximization approach to blind separation and blind deconvolution[END_REF]. The equivalence between the InfoMax principle and ML was proven in [Cardoso, 1997, Obradovic and[START_REF] Obradovic | [END_REF]]. Many algorithms were proposed thereafter. Pham and Garat provided two methods by using prior distributions of the sources [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF]. Boscolo et al. proposed an ICA method based on a nonparametric kernel density estimation technique [START_REF] Boscolo | Independent component analysis based on nonparametric density estimation[END_REF]. Dégerine and Zaïdi derived an ML approach for Gaussian autoregressive sources [START_REF] Dégerine | Separation of an instantaneous mixture of Gaussian autoregressive sources by the exact maximum likelihood approach[END_REF]. Kokkinakis and Nandi proposed to use generalized gamma densities in the ML estimation for ICA [START_REF] Kokkinakis | Generalized gamma density-based score functions for fast and flexible ICA[END_REF]. Xue et al. solved the ML-based ICA by means of gradient equation [START_REF] Xue | Independent component analysis based on gradient equation and kernel density estimation[END_REF]. Ge and Ma studied the spurious solutions of the ML approach [START_REF] Ge | Spurious solution of the maximum likelihood approach to ICA[END_REF]. Vía et al. provided an ICA method for quaternion signals [START_REF] Vía | Quaternion ICA from second-order statistics[END_REF].

ICA by minimization of mutual information

Mutual Information (MI) is a well-established information-theoretic measure of statistical dependence, based on which ICA approaches can be developed. MI between P random variables, s i , i ∈ {1, . . . , P }, is defined as follows:

I(s 1 , s 2 , . . . , s P ) = P i=1 H(s i ) -H(s) (1.7)
where the negentropy function H(•) is defined in equation (1.4). MI is nonnegative and equal to zero if and only if the random variables are statistically independent. MI-based ICA methods intend to minimize the MI of the separated sources s = V x with respect to the demixing matrix V , yielding the maximum independence. MI can serve as a unifying framework for maximization of non-Gaussianity and ML estimation [Yang, 1997, Obradovic and[START_REF] Obradovic | [END_REF]. A large number of ICA methods are based on the minimization of MI. He et al. proposed a nonparametric algorithm by means of clustering-based multivariate density estimation [START_REF] He | Blind source separation using clustering-based multivariate density estimation algorithm[END_REF]. Pham proposed a method based on higher order statistics [Pham, 2000], and a method for colored sources via the Gaussian MI criterion [Pham, 2001a]. Babaie-Zadeh et al. designed an algorithm based on differential MI [START_REF] Babaie-Zadeh | Differential of the mutual information[END_REF]. Kraskov et al. gave several MI estimators for ICA [START_REF] Kraskov | Estimating mutual information[END_REF]. Babaie-Zadeh and Jutten introduced a nonparametric gradient method [START_REF] Babaie-Zadeh | A general approach for mutual information minimization and its application to blind source separation[END_REF]. A fast fixed-point algorithm was provided by Hulle [Hulle, 2008]. Other methods include [Almeida, 2003,Pham, 2004, Li and Adalı, 2010, Rashid and Yu, 2012], to cite a few.

Independent component analysis

Chapter 1

ICA by tensor decomposition of cumulant arrays

Higher-Order (HO) tensors are essentially HO arrays whose components have more than two indices, and they can be seen as HO extensions of vectors (first order arrays) and matrices (second order arrays). HO cumulant arrays can be considered as HO generalizations of the covariance matrix, and they are the measure of statistical dependence at higher order. An important class of ICA approaches, as well as the proposed methods in this thesis report, consists of using HO cumulant arrays.

Preliminaries

Now we introduce some basic definitions in multilinear algebra which are necessary for the cumulant-based ICA methods.

Definition 1.1. The outer product

T = u (1) • u (2) • • • • • u (q) of q vectors u (i) ∈ R N i , i ∈ {1, . . . , q}, is a q-th order array in R N 1 ×N 2 ו••×Nq whose entries are defined by T i 1 ,i 2 ,...,iq = u (1) i 1 u (2) i 2 • • • u (q) iq .
Definition 1.2. Each q-th order array T expressed as the outer product of q vectors is a rank-1 q-th order array.

More generally, the rank of a q-th order array is defined as follows:

Definition 1.3. The rank of an array T ∈ R N 1 ×N 2 ו••×Nq , denoted by rank(T ), is the minimal number of rank-1 arrays belonging to R N 1 ×N 2 ו••×Nq that yield T in a linear combination.
Despite the similarity between the definition of the tensor rank and its matrix counterpart, the rank of a HO array may exceed its dimensions [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]. An important concept of matrix rank named Kruskal rank is defined as follows [START_REF] Harshman | The PARAFAC model for three-way factor analysis and multidimensional scaling[END_REF]:

Definition 1.4. The Kruskal rank of a matrix U , denoted by rank K (U ), is defined as the largest number such that every subset of columns of U is linearly independent.

We can extract a matrix slice from a HO array by means of the following definition: Definition 1.5. A q-th order array matrix slice is a 2-dimensional section (fragment) of a q-th order array, obtained by fixing q -2 of the q indices [Cichocki et al., 2009]. Now we define the multiplication of a HO array with a matrix: Definition 1.6. The mode-i product of an array

T ∈ R N 1 ×N 2 ו••×Nq by a matrix U ∈ R J i ×N i , denoted by T × i U is an array in R N 1 ×N 2 ו••×N i-1 ×J i ×N i+1 ו••×Nq defined by: (T × i U ) n 1 ,n 2 ,...,j i ,...,nq = N i n i =1 T n 1 ,n 2 ,...,n i ,...,nq U j i ,n i (1.8)
for all index values of j i .

According to definitions 1.1, 1.2 and 1.3, the low-rank Canonical Polyadic Decomposition (CPD) model of a given HO array is defined by: Chapter 1 1.2. Independent component analysis Definition 1.7. Given an integer number P of components, the CPD model of a q-th order array T ∈ R N 1 ×N 2 ו••×Nq can be expressed as follows:

T = P p=1 u (1) p • u (2) p • • • • • u (q) p + R (1.9)
where the q-th order array R represents the model residual.

The notation T = [[U (1) , U (2) , . . . , U (q) ]] + R refers to the CPD of T by equation (1.9) with the associated loading matrices

U (i) = [u (i) 1 , • • • , u (i) P ] ∈ R N i ×P . When R is a null
array, we have an exact CPD. An exact CPD is considered to be essentially unique when it is only subject to scale and permutation indeterminacies. A sufficient condition of uniqueness was defined by means of Kruskal rank as presented in the following theorem [START_REF] Sidiropoulos | On the uniqueness of multilinear decomposition of n-way arrays[END_REF]]:

Theorem 1.1. The CPD of a rank-P q-th order array T = [[U (1) , U (2) , . . ., U (q) ]] is considered to be essentially unique when the following inequality holds:

q i=1 rank K (U (i) ) ≥ 2P + (q -1) (1.10)
More refined uniqueness conditions specifically for third-order arrays were addressed in [Kruskal, 1977, Ten Berge and Sidiropoulos, 2002, Jiang and Sidiropoulos, 2004, De Lathauwer, 2006].

The CPD problem and its algorithms have been studied extensively for the past four decades, see [Harshman, 1970, Harshman and Lundy, 1994, Comon, 2002, De Lathauwer et al., 2004, De Lathauwer, 2006, Rajih et al., 2008, Comon et al., 2009, Kolda and Bader, 2009, Karfoul et al., 2011, Röemer and Haardt, 2013, Hajipour Sardouie et al., 2013, Luciani and Albera, 2014, Comon, 2014] and references therein.

As a special form of the CPD of a third-order array where two loading matrices are equal, the low-rank INdividuals Differences in SCALing (INDSCAL) analysis [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of Eckart-Young decomposition[END_REF] is defined as follows:

Definition 1.8. For a given P , corresponding to the number of rank-1 terms, the IND-SCAL decomposition of a third-order array T ∈ R N 1 ×N 1 ×N 3 can be expressed as:

T = P p=1 u (1) p • u (1) p • u (3) p + R (1.11)
with the loading matrices

U (1) = [u (1) 1 , • • • , u (1)
P ] ∈ R N 1 ×P and U (3) = [u (3) 1 , • • • , u (3)
P ] ∈ R N 3 ×P , and R being the residual term.

The INDSCAL decomposition can also be described by using the frontal slices of T :

∀ k ∈ {1, 2, • • • , N 3 }, T (k) = T :,:,k = U (1) D (k) U 3 (U (1) ) T + R (k) (1.12)
where T (k) = T :,:,k is the k-th frontal slice of T , where D

U 3 ∈ R P ×P is a diagonal matrix whose diagonal contains the elements of the k-th row vector of U (3) , and where R (k) = R :,:,k is the k-th frontal slice of R. In fact, equation (1.12) shows that the 1.2. Independent component analysis

Chapter 1

INDSCAL decomposition problem is essentially a joint diagonalization by congruence problem which we will address in the next section. The uniqueness condition for the INDSCAL decomposition was established by Afsari [Afsari, 2008]. The q-th order cumulant array of a real random vector is essentially a q-th order array, and it can be related to moments of order lower than and equal to q using Leonov-Shiryaev formula [McCullagh, 1987]:

Property 1.1. The q-th order cumulant array C (q)

x ∈ R N ×N ו••×N of a real random vector x ∈ R N satisfies the following element-wise equation:

C (q) x,n 1 ,n 2 ,...,nq = cum{x n 1 , x n 1 , . . . , x nq } = (-1) -1 ( -1)! E{ n∈I 1 x n } E{ n∈I 2 x n } • • • E{ n∈I x n } (1.13)
where the summation includes all possible partitions {I 1 , I 2 , . . . , I } (1 ≤ ≤ q) of the integers {n 1 , n 2 , . . . , n q } (see [START_REF] Karfoul | Blind underdetermined mixture identification by joint canonical decomposition of HO cumulants[END_REF] for more details).

For a real zero-mean random vector x, the cumulants up to order four are expressed in details as follows:

C (1) x,n 1 = E{x n 1 } (1.14) C (2) x,n 1 ,n 2 = E{x n 1 x n 2 } (1.15) C (3) x,n 1 ,n 2 ,n 3 = E{x n 1 x n 2 x n 3 } (1.16) C (4) x,n 1 ,n 2 ,n 3 ,n 4 = E{x n 1 x n 2 x n 3 x n 4 } -E{x n 1 x n 2 } E{x n 3 x n 4 } -E{x n 1 x n 3 } E{x n 2 x n 4 } -E{x n 1 x n 4 } E{x n 2 x n 3 } (1.17)
HO cumulant arrays enjoy a number of important properties, which are listed as follows [START_REF] Nikias | [END_REF]Mendel, 1993, De Lathauwer, 2010]:

1. Symmetry: real cumulant arrays are fully symmetric under an arbitrary permutation of the indices.

2. Multi-linearity: if a real random vector x admits the linear transformation x = As, then the following equality holds:

C (q) x = C (q) s × 1 A × 2 A • • • × q A (1.18)
3. Partitioning of independent variables: if a subset of p random variables x n 1 , x n 2 , . . ., x np is statistically independent of the other variables, then their HO cumulants are equal to zero. Therefore, a HO cumulant array of a random vector that has mutually independent components is a diagonal array. In this case only the entries of which all the indices are equal can be different from zero.

4. Sum of independent variables: the cumulant array of a sum of independent random vectors is the sum of the individual cumulant arrays.

5. Non-Gaussianity: HO cumulant arrays of a Gaussian variable are zero arrays. In combination with the fourth property, for a random vector x corrupted by an additive Gaussian noise vector ε, we have, for q > 2:

C (q) x+ε = C (q) x + C (q) ε = C (q) x (1.19)
Chapter 1

1.2. Independent component analysis Consequently, we can state that HO cumulants of a set of random variables give an indication of their mutual statistical dependence, and that the HO cumulant of a single random variable is some measure of its non-Gaussianity [START_REF] Pesquet | [END_REF]Moreau, 2001, De Lathauwer, 2010]. Moreover, the above properties of a HO cumulant array determine the algebraic structure of that array, which establishes the algorithmic foundation of a family of well-known ICA methods and also of the proposed algorithms in this thesis report.

Cumulant-based ICA algorithms

Considering the ICA model (1.2) and using the properties of cumulants discussed in the previous section, the q-th order cumulant array of the observation vector x can be expressed as follows:

C (q) x = C (q) s × 1 A × 2 A • • • × q A + C (q) ε = P p=1 C (q) s,p,...,p a p • a p • • • • • a p + C (q) ε (1.20)
where a p is the p-th column vector of the mixing matrix A, where C

(q) ε is the q-th order cumulant array of the noise vector, and where C (q) s,p,...,p denotes the q-th order cumulant of the p-th source. The second equality holds when the sources are statistically independent, in which case the cumulant array C (q) s of the source vector is diagonal. Moreover, C (q) ε (q > 2) vanishes if the noise has a Gaussian distribution. Equation (1.20) is essentially the CPD of a cumulant array. Under the uniqueness condition of the CPD on the matrix A, and assuming that all the sources have non-zero q-th order cumulants C (q) s,p,...,p and ε is a Gaussian noise vector, the mixing matrix A can be estimated uniquely by CPD of C (q)

x (q > 2) up to scale and permutation indeterminacies. Cardoso proposed a super-symmetric decomposition method for the Fourth-Order (FO) cumulant array [Cardoso, 1991]. Comon provided a CPD method for cumulant arrays but restricted the mixing matrix A to be of size (2 × 3) [Comon, 2004]. Albera et al. proposed an algorithm to compute the CPD of a Hermitian FO array and applied it the FO cumulant array of the observation array in order to identify A [START_REF] Albera | ICAR: a tool for blind source separation using fourth-order statistics only[END_REF]. The method was extended to 2q-th order Hermitian arrays by the same authors in [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF], allowing for the identification of underdetermined mixing matrices (i.e. with more columns than rows). Ferréol et al. showed how to combine several time-delayed FO cumulant arrays for the identification of underdetermined mixing matrices, provided that one of the FO cumulant arrays was positive (or negative) semi-definite [START_REF] Ferréol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF]. De Lathauwer and Castaing proposed to decompose jointly several Second-Order (SO) time-delayed cumulant arrays [De Lathauwer and Castaing, 2008]. Karfoul et al. provided a joint CPD method of several HO cumulant arrays relaxing the positive (or negative) semi-definite assumption [START_REF] Karfoul | Blind underdetermined mixture identification by joint canonical decomposition of HO cumulants[END_REF]. A generalization of the Enhanced Line Search Alternating Least Squares (ELS-ALS) method [START_REF] Rajih | Enhanced line search: A novel method to accelerate PARAFAC[END_REF] to the case of complex-valued arrays of order strictly greater than 3 and a new algorithm for the CPD of Hermitian HO arrays were proposed in [START_REF] Karfoul | Iterative methods for the canonical decomposition of multi-way arrays: Application to blind underdetermined mixture identification[END_REF].

On the other side, considering the cumulant array of the estimated source vector s = V x: s into a diagonal array, yielding s having statistically independent components. Suppose that the cumulant array is computed on the whitened observation vector. The search space of V is then reduced to the group of (P × P ) unitary matrices. A (P × P ) unitary matrix can be parameterized as a product of P (P -1)/2 Givens rotation matrices. Comon proposed to use the FO cumulant array and to estimate these Givens rotation matrices sequentially with Jacobi iterations. In each Jacobi iteration, the parameter of one rotation matrix can be estimated either by maximizing the total sum of squares of the diagonal components of the cumulant array C (q) s [Comon, 1994], or by maximizing the trace of that cumulant array if all the source cumulants are known to be positive [START_REF] Comon | Improved contrast dedicated to blind separation in communications[END_REF]. The former algorithm is known as CoM2, and the latter is called CoM1, where the numbers "2" and "1" indicate that the algorithms maximize the 2-norm and 1-norm of the diagonal components, respectively. Cardoso and Souloumiac proposed to maximize another cost function (called contrast [Comon, 1994, Moreau, 2001]) based on FO cumulants. This maximization problem was reformulated as the Joint Approximate Diagonalization of Eigen-matrices, leading to the JADE algorithm [START_REF] Cardoso | Blind beamforming for non-Gaussian signals[END_REF]. On the other hand, De Lathauwer et al. considered the FO cumulant array as a stack of N Third Order (TO) arrays of dimension (N × N × N ), and proposed an ICA algorithm based on Simultaneous Third-Order Tensor Diagonalization (STOTD) [De Lathauwer et al., 2001]. Moreau extended the JADE and STOTD methods to any HO cumulant by providing a new family of contrast functions that served as a unifying framework for JADE and STOTD [Moreau, 2001]. Moreover, he revealed a link between these new contrast functions and a joint diagonalization criterion, yielding a novel method called eJADE. Tobias and Wiskott extended the CoM2 algorithm by incorporating both the TO and FO cumulant arrays [START_REF] Blaschke | CuBICA: independent component analysis by simultaneous third-and fourth-order cumulant diagonalization[END_REF]]. However, it was shown in [Moreau, 2006] that the contrast function used by Tobias and Wiskott was a particular case of the general contrast introduced by Moreau in [Moreau, 2001]. Another ICA method called Second Order Blind Identification (SOBI) for separating stationary sources with different spectral contents was proposed by Belouchrani et al., which aimed at jointly diagonalizing a set of whitened SO time-delayed cumulant arrays [Belouchrani et al., 1997].

C (q) s = C (q) x × 1 V × 2 V • • • × q V (1.
As informed in the previous paragraphs, Joint Diagonalization by Congruence (JDC) is the core step in quite a number of cumulant-based ICA methods, such as JADE [START_REF] Cardoso | Blind beamforming for non-Gaussian signals[END_REF], eJADE [Moreau, 2001] and SOBI [Belouchrani et al., 1997]. The connection between ICA and JDC, as well as the JDC algorithms, will be addressed in the next section.

From ICA to joint diagonalization by congruence

Joint diagonalization of cumulant matrix slices by congruence is a principal procedure for many ICA algorithms [START_REF] Moreau | Blind Identification and Separation of Complex-Valued Signals[END_REF]. In this section, firstly we study the matrix-wise structure of cumulant arrays and show that the matrix slices of cumulant arrays preserve a joint congruence transformation structure. The mixing matrix of ICA can be estimated by means of a JDC technique. Secondly, the uniqueness and identifiability issues are briefly addressed. Thirdly, we provide an overview of existing JDC algorithms.

Chapter 1 1.3. From ICA to joint diagonalization by congruence

JDC problem

A matrix slice of a q-th order cumulant array C (q)

x (q > 2) can be obtained by varying its first two indices and picking one particular value for the other indices. If x obeys the ICA model (1.2) and is contaminated by a Gaussian noise, from equation (1.20), the matrix slices have the following structure:

C (q)
x,:,:,n 3 ,...,nq =

P n 1 =1 (C (q) s,n 1 ,...,n 1 A n 3 ,n 1 • • • A nq,n 1 ) a n 1 • a n 1 = A D (n 3 ,...,nq) A T (1.22)
where A n i ,n j is the (n i , n j )-th component of A, where a n 1 denotes the n 1 -th column vector of A, where D (n 3 ,...,nq) ∈ R P ×P are diagonal matrices, and where the equality .22) shows that the matrix slices of cumulant arrays share the following joint congruence transformation: K} (1.23) where D (k) ∈ R P ×P are diagonal matrices, and where K = N q-2 is the total number of matrix slices from a q-th order cumulant array. In this case, K matrices C (k) ∈ R N ×N are jointly diagonalizable and A ∈ R N ×P is called the joint transformation matrix. It is noteworthy that in order to use HO cumulant arrays and their matrix slices in ICA, the sources are assumed to be mutually independent and non-Gaussian. Consequently, in order to identify the mixing matrix A from the cumulant matrix slices, we aim at solving the following JDC problem: Problem 1.2. Given K symmetric matrices C (k) ∈ R N ×N and an integer P , find a joint transformation matrix A ∈ R N ×P and K diagonal matrices D (k) ∈ R P ×P such that these matrices satisfy equation (1.23).

a n 1 • a n 1 = a n 1 a T n 1 is used implicitly. Equation (1
C (k) = AD (k) A T , k ∈ {1, . . . ,
In most papers, the terms "joint diagonalization" and "simultaneous diagonalization" were utilized to state the JDC problem, such as in [Yeredor, 2002]. We propose to emphasis the "Congruence" transformation in the Joint Diagonalization problem, leading to the acronym "JDC", in oder to avoid confusion with the joint diagonalization by equivalence problem C (k) = AD (k) A -1 , which is an essentially different problem when A is not unitary [START_REF] Luciani | Joint eigenvalue decomposition using polar matrix factorization[END_REF],Luciani and Albera, 2011,Luciani and Albera, 2014].

Uniqueness and identifiability

If the ICA model (1.2) is noiseless and the mixing matrix A ∈ R N ×P is of full column rank with N > P , one can always find P of the N observations which are linearly independent and discard the other observed mixtures without loss of information, leading to a square mixing matrix. More generally, an appropriate (P × N ) matrix transformation can be applied to the N observations in order to merge them and to reduce the dimension of the observation vector. Therefore, without loss of generality, we can assume that for a noiseless JDC problem the joint transformation matrix A is a nonsingular square matrix of dimension (N × N ).

Afsari proposed a uniqueness condition for the noiseless JDC problem by measuring the coherence of the K diagonal matrices D (k) [Afsari, 2008]. Let D be a (K ×N ) matrix 1.3. From ICA to joint diagonalization by congruence Chapter 1 whose (k, n)-th component is defined as follows:

D k,n = D (k) n,n , k ∈ {1, . . . , K}, n ∈ {1, . . . , N } (1.24)
where

D (k)
n,n is the (n, n)-th component of the k-th diagonal matrix D (k) . Let d n 1 and d n 2 denote the n 1 -th and n 2 -th column vectors of D, respectively. The cosine of angle ψ n 1 ,n 2 between two vectors d n 1 and d n 2 can be derived by using the following Euclidean dot product formula:

d T n 1 d n 2 = d n 1 d n 2 cos(ψ n 1 ,n 2 ) (1.25)
Definition 1.9. The coherence ρ of K diagonal matrices D (k) is defined as the maximum absolute cosine of angle ψ n 1 ,n 2 between the column vectors of D as follows:

ρ = max n 1 ,n 2 n 1 =n 2 | cos(ψ n 1 ,n 2 )| with cos(ψ n 1 ,n 2 ) = d T n 1 d n 2 d n 1 d n 2 (1.26)
If one of the two vectors is zero, ρ is defined to be one by convention.

The coherence ρ, also known as the modulus of uniqueness, captures the uniqueness property of the JDC problem as described in the following theorem [Afsari, 2008].

Theorem 1.2. Let each C (k) satisfy equation (1.23). The necessary and sufficient condition to identify a unique nonorthogonal joint transformation matrix A up to scale and permutation indeterminacies is ρ < 1.

The uniqueness issue of the complex-valued JDC problem was addressed in [START_REF] Kleinsteuber | Uniqueness analysis of non-unitary matrix joint diagonalization[END_REF]. Based on theorem 1.2 and the properties of cumulants, Afsari provided an identifiability condition for the JDC-based ICA method [Afsari, 2008].

Theorem 1.3. (Identifiability of ICA -JDC formulation) Regarding the noiseless ICA model, let us assume that i) the covariance matrix R s of the source vector s is nonsingular, ii) the q-th order marginal cumulants (for some q > 2) of all the source components exist and iii) at most one of them is zero. Then an exact JDC of the N q-2 cumulant matrix slices C (q)

x,:,:,n 3 ,...,nq as well as the covariance matrix R x , results in finding A up to column scaling and permutation. For a source vector with finite cumulants of all orders, this process fails to identify A if and only if the source vector has more than one Gaussian component.

Note that the use of the covariance matrix is optional when all the source components have non-zero q-th order cumulants.

JDC algorithms

The JDC problem has been widely studied during the past two decades and has been mostly handled as an optimization problem. The algorithms mainly depend on the criterion chosen to perform the optimization. These criteria include the indirect-fit criterion, the direct-fit criterion, the combination of the direct and indirect-fit criteria and the approximate log-likelihood criterion. Following the recent survey of Chabriel et al. [START_REF] Chabriel | Joint matrices decompositions and blind source separation: A survey of methods, identification, and applications[END_REF], we provide a comprehensive overview of these four groups of criteria and the associated JDC algorithms. In the rest of this section, when a JDC problem is mentioned without further specification, it refers to a "square" problem, in which the matrices A, D (k) and C (k) are of size (N × N ).

Chapter 1 1.3. From ICA to joint diagonalization by congruence

Minimizing the indirect-fit criterion

Based on the JDC model (1.23), the matrices D (k) can be expressed as a function of A by D (k) = A -1 C (k) (A -1 ) T and should be diagonal by definition. A natural way to compute A is to minimize the squares of the off-diagonal components of D (k) with respect to A, leading to the following indirect-fit criterion:

Ψ 1 (A) = K k=1 off(A -1 C (k) A -T ) 2 F (1.27)
where off(•) sets the diagonal entries of the matrix argument to zero, and where A -T denotes (A -1 ) T . Frequently the minimization of criterion (1.27) is performed with respect to a matrix V such that V = A -1 , instead of A for simplicity, and V is called the joint diagonalizer. In the context of ICA, V can be regarded as a demixing matrix. To use this criterion the matrix V (or A -1 ) should be properly constrained in order to avoid the trivial zero solution and/or degenerate solutions. In the following, we present the options of constraints and penalty terms appeared in literature:

1. V is unitary. Cardoso and Souloumiac proposed an efficient algorithm, namely Joint Approximate Diagonalization (JAD), which estimated V by successive Jacobi iterations [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF]. Moreau extended this algorithm by introducing a generalized JDC criterion [Moreau, 2001]. Zhang et al. estimated V by sequential Householder transforms [START_REF] Zhang | Joint approximate diagonalization using bilateral rank-reducing Householder transform with application in blind source separation[END_REF].

2. Each column vector of V has unit norm. This constraint is less restrictive than the former one and was adopted in [START_REF] Wang | A new joint diagonalization algorithm with application in blind source separation[END_REF],Dégerine and Kane, 2007,Fadaili et al., 2007].

3. The diagonal elements of the matrix V C (k 0 ) V T must be equal to one, where C (k 0 ) with k 0 ∈ {1, . . . , K} is a positive definite matrix in {C (k) }. In ICA, C (k 0 ) could be the zero-delayed covariance matrix of the observations, and this constraint turns out to force the separated sources having unit variance. Vollgraf and Obermayer proposed a Quadratic DIAGonalization algorithm (QDIAG), which estimated each row vector of V alternately by matrix eigenvector computation [START_REF] Vollgraf | Quadratic optimization for simultaneous matrix diagonalization[END_REF].

4. V has a unit determinant. This constraint can be guaranteed by expressing V as a product of a sequence of elementary unit-determinant matrices, such as Givens rotation matrices, hyperbolic rotation matrices and elementary triangular matrices. A merit of using these elementary matrices is that each matrix is determined by only one parameter, which generally can be estimated algebraically by means of a Jacobi-like optimization procedure. These methods include the LUJ1D and QRJ1D algorithms proposed by Afsari [Afsari, 2006], the J-DI algorithm proposed by Souloumiac [Souloumiac, 2009], the SL algorithm proposed by Sørensen et al. [START_REF] Sørensen | Approximate tensor diagonalization by invertible transforms[END_REF], and the coupled LU and QR factorizations based algorithms proposed by Maurandi et al. [Maurandi et al., 2013] and by [START_REF] Maurandi | A new fast Jacobi-like algorithm for non-orthogonal joint diagonalization of real-valued matrices based on a QR parameterization[END_REF], respectively. Several methods extended the ideas for complex-valued JDC problems [START_REF] Guo | Approximate joint diagonalization by nonorthogonal nonparametric Jacobi transformations[END_REF], Gong et al., 2012, Mesloub et al., 2014, Maurandi and Moreau, 2014a, Maurandi et al., 2014], to cite a few.

1.3. From ICA to joint diagonalization by congruence Chapter 1

5. V is constrained to be invertible. Ziehe et al. proposed an iterative method, namely Fast Frobenius DIAGonalization (FFDIAG) [Ziehe et al., 2004]. Given a nonsingular initialization V (0) , in the it-th iteration, V (it-1) is updated by a strictly diagonally dominant matrix U (i.e. for every row of U , the magnitude of the diagonal entry in a row is larger than the sum of the magnitudes of all the other non-diagonal entries in that row), therefore the product V (it) = U V (it-1) is guaranteed to be invertible. A complex-valued extension of FFDIAG was presented in [START_REF] Xu | A fast algorithm for nonunitary joint diagonalization and its application to blind source separation[END_REF]. Yeredor et al. derived a relative gradient approach and pointed out that for small enough step-sizes the invertibility of V was ensured [START_REF] Yeredor | Approximate joint diagonalization using a natural gradient approach[END_REF]. Thereafter, Ghennioui et al. improved the former algorithm based on computing an optimal step-size [START_REF] Ghennioui | An optimal step size relative gradient based joint diagonalization algorithm[END_REF]. Trainini et al. extended the relative gradient algorithm for complex-valued matrices [START_REF] Trainini | A relative gradient algorithm for joint decompositions of complex matrices[END_REF].

Several variations of gradient-like methods were provided by [START_REF] Trainini | Variations around gradient like algorithms for joint diagonalization of hermitian matrices[END_REF].

6. Several methods resorted to add penalty terms to criterion (1.27). In [START_REF] Wang | Penalty function-based joint diagonalization approach for convolutive blind separation of nonstationary sources[END_REF], Wang et al. introduced two penalty terms, namely

V V T -I N 2 F and diag(V -I N ) 2
F , and proposed a gradient descent method, where diag(.) returns a matrix comprising only the diagonal elements of the input matrix. Li and Zhang utilized log | det(V )| and performed the minimization with respect to each column vector of V alternately [Li and Zhang, 2007]. For a non-square matrix V , Li et al. turned to use log | det(V V T )| and proposed a gradient decent algorithm [Li et al., 2007]. This penalty term was also adopted by [START_REF] Zhou | Nonorthogonal approximate joint diagonalization with well-conditioned diagonalizers[END_REF]. Pham and Congedo proposed a new intrinsic constraint k (V n,: C (k) V T n,: ) 2 = 1, and derived a pseudo Newton algorithm, where V n,: is the n-th row vector of V [START_REF] Pham | Least square joint diagonalization of matrices under an intrinsic scale constraint[END_REF].

Besides the criterion Ψ 1 (A) in (1.27), Afsari proposed another indirect-fit criterion suitable for JDC, which was scale-invariant in A (and V ) [Afsari, 2006, Afsari, 2008]:

Ψ 2 (A) = K k=1 C (k) -A diag(A -1 C (k) A -T )A T 2 F (1.28)
The scale invariance means that the change of the scale of any columns of A does not affect the value of this criterion. Based on criterion (1.28), Afsari proposed two algorithms, namely LUJ2D and QRJ2D, by means of LU and QR matrix factorizations of V , respectively [Afsari, 2006]. Recently, Yeredor et al. modified this criterion and used it in a different JDC problem, named sequentially drilled JDC [START_REF] Yeredor | A sequentially drilled joint congruence (SeDJoCo) transformation with applications in blind source separation and multiuser MIMO systems[END_REF].

Minimizing the direct-fit criterion

The direct-fit criterion is a measure of the squared difference between the given matrices C (k) and their assumed model parameters in terms of estimating the joint transformation matrix A and K diagonal matrices

D (k) (k ∈ {1, . . . , K}), Ψ 3 (A, {D (k) }) = K k=1 C (k) -AD (k) A T 2 F (1.29)
Chapter 1 1.3. From ICA to joint diagonalization by congruence This criterion was originally introduced by Wax and Sheinvald [START_REF] Wax | A least-squares approach to joint diagonalization[END_REF]. They revealed that criterion (1.29) coincided with criterion (1.27) when A was constrained to be unitary. For a non-unitary matrix A, Van der Veen proposed to minimize (1.29) via subspace fitting technique [ Van der Veen, 2001]. An effective optimization scheme was proposed by Yeredor [Yeredor, 2002]. The minimization of (1.29) was performed alternately between each column vector of A and the diagonal matrices D (k) , giving birth to the well-known Alternating Columns and Diagonal Center (ACDC) algorithm. Recently, Zhang and Lou developed a recursive algorithm, which provided a compromise between performance and numerical complexity [START_REF] Zhang | A recursive solution to nonunitary joint diagonalization[END_REF]. Zeng and Feng proposed a novel algorithm based on the hybrid trust region method [START_REF] Zeng | Non-orthogonal joint diagonalization algorithm based on hybrid trust region method and its application to blind source separation[END_REF]. In addition, criterion (1.29) does not require A to be square, and it can even handle the underdetermined case where A has more columns than rows [De Lathauwer and Castaing, 2008].

On the other hand, in the absence of noise, A can be extracted from the Generalized EigenValue Decomposition (GEVD) of two matrices C (1) and C (2) as follows:

C (1) (C (2) ) -1 = A(D (1) (D (2) ) -1 )A -1 (1.30)
where we assume that D (2) (and hence C (2) ) is nonsingular. When a noisy JDC problem involving more than two matrices is considered, Yeredor provided a suboptimal but closed-form (non-iterative) solution by deriving two "representative" matrices from {C (k) } and computing the GEVD [Yeredor, 2005]. The algorithm is called Exact Joint Diagonalization (EJD). Based on EJD, Chabriel and Barrère proposed a more robust method, namely DIagonalization using Equivalent Matrices (DIEM) [START_REF] Chabriel | A direct algorithm for nonorthogonal approximate joint diagonalization[END_REF]. DIEM requires that at least K = N matrices C (k) are available.

Additionally, each matrix C (k) ∈ R N ×N can be regarded as a fontal matrix slice of a TO array C ∈ R N ×N ×K . The minimization of criterion (1.29) is directly linked

to the CPD of C, such that C = [[A, A, D]],
where D is defined in equation (1.24). As previously reviewed in section 1.2.3.1, this special case of CPD when two loading matrices coincide is called INDSCAL decomposition [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of Eckart-Young decomposition[END_REF]. Comon and Rajih proposed a method for INDSCAL decomposition in [Comon and Rajih, 2006, Appendix 7.2]. Albera and Karfoul derived a different Alternating Least Squares (ALS) algorithm (ALS) incorporating an Enhanced Line Search (ELS) procedure [START_REF] Albera | A fast INDSCAL analysis for blind underdetermined mixture identification from several even higher order cumulants[END_REF]. Trainini and Moreau proposed another ELS-ALS method for complexvalued matrices [START_REF] Trainini | A least squares algorithm for global joint decomposition of complex matrix sets[END_REF]. A damped Gaussian-Newton algorithm was provided by [START_REF] Koldovský | Stability analysis and fast damped-Gauss-Newton algorithm for INDSCAL tensor decomposition[END_REF]. Moreover, since INDSCAL model is a special case of the CPD, one can attempt to decompose C by general CPD C = [[A (1) , A (2) , D]], and suppose that eventually the two loading matrices A (1) and A (2) will converge to an estimate of A. But in practice, some differences will remain between A (1) and A (2) . Then the n-th column vector of A, denoted by a n , can be estimated as the dominant left singular vector of the matrix [a

(1) n a (2) n ] ∈ R N ×2 , where a (1)
n and a

(2) n are the n-th column vectors of A (1) and A (2) , respectively, and n ∈ {1, . . . , N } [De [START_REF] Lathauwer | [END_REF]. An explicit link between the CP and JDC decompositions was established by De Lathauwer in [De Lathauwer, 2006].

From ICA to joint diagonalization by congruence

Chapter 1

Combination of the direct and indirect-fit criteria

In criterion (1.29), instead of fitting the set of matrices C (k) by A and D (k) directly, we can consider fitting a set of partially diagonalized matrices V (it) C (k) (V (it) ) T at the it-th iteration by means of a residual joint transformation Ā(it) and K diagonal matrices D (k) v , using the following criterion:

Ψ 4 ( Ā(it) ) = K k=1 V (it) C (k) (V (it) ) T -Ā(it) D (k) v ( Ā(it) ) T 2 F (1.31)
where it denotes the iteration index, where V (it) is the obtained joint diagonalizer at the it-th iteration, and where

D (k) v = diag{V (it) C (k) (V (it) ) T }.
Criterion (1.31) can be regarded as the combination of the direct criterion (1.27) and the indirect-fit criterion (1.29). Each minimization step consists in one iteration of a direct-fit procedure with respect to Ā(it) . In the direct-fit procedure, a Gaussian-Newton method is adopted for the purpose of fast convergence in a neighborhood of the true local minimum which is close to Ā(it) = I N . By setting I N as the initial point, only one iteration of the Gaussian-Newton optimization is applied for each step, because at such an initial point, the inverse of the Hessian matrix can be computed by solving distinct sets of (2 × 2) linear equations [START_REF] Chabriel | Joint matrices decompositions and blind source separation: A survey of methods, identification, and applications[END_REF]. Once the residual joint transformation matrix it) . This method was proposed by Tichavský and Yeredor, and called Uniformly Weighted Exhaustive Diagonalization with Gauss itErations (UWEDGE) [START_REF] Tichavský | Fast approximate joint diagonalization incorporating weight matrices[END_REF].

Ā(it) is found, the joint diagonalizer V (it+1) is computed by V (it+1) = ( Ā(it) ) -1 V (

Minimizing the approximate log-likelihood criterion

As already discussed in section 1.3.1, when the ICA model assumes statistically independent sources with variances that are constant within each block but varying between blocks, the jointly diagonalizable matrices C (k) can be sample covariance matrices taken from distinct blocks of the observations [START_REF] Pham | Blind separation of instantaneous mixtures of nonstationary sources[END_REF]. Consequently, each C (k) is a positive definite matrix. Pham utilized the following criterion, which measured the diagonality of positive definite matrices [START_REF] Pham | Joint approximate diagonalization of positive definite Hermitian matrices[END_REF]:

Ψ 5 (A) = K k=1 log det {diag(A -1 C (k) A -T )} det {A -1 C (k) A -T } (1.32)
Minimizing (1.32) is approximately equivalent to maximizing a log-likelihood objective function for JDC. Historically, this criterion was first studied by Flury and Gautsch for seeking a unitary joint transformation matrix A [START_REF] Flury | An algorithm for simultaneous orthogonal transformation of several positive definite symmetric matrices to nearly diagonal form[END_REF]. In this case, the denominators in (1.32) are constant. Flury and Gautsch estimated A by successive Jacobi iterations. Pham considered the more general case that A was non-unitary.

He proposed to maximize a lower bound of the decrease of (1.32) by means of sequential (2 × 2) transformations, and derived algebraic solutions for those transformation matrices [START_REF] Pham | Joint approximate diagonalization of positive definite Hermitian matrices[END_REF]. Thereafter, Joho provided a Newton algorithm [Joho, 2008].

Todros and Tabrikian proposed a computationally efficient algorithm by minimizing an approximation of criterion (1.32) [START_REF] Todros | QML-based joint diagonalization of positive-definite Hermitian matrices[END_REF].

Chapter 1

Nonnegative matrix factorization

In many applications, the performance of the JDC algorithms can be substantially improved by introducing appropriate weighting in the optimization criterion. For example, the indirect-fit criterion (1.27) can be rewritten equivalently as follows:

Ψ 1 (A) = [f(A)] T f(A) (1.33)
where the function f(A) returns a KN (N -1)/2-dimensional column vector composed of all off-diagonal components below the main diagonal of the

K matrices A -1 C (k) A -T .
The weighted indirect-fit criterion can be derived by inserting a positive definite weight matrix Υ in (1.33) as follows [START_REF] Tichavský | Fast approximate joint diagonalization incorporating weight matrices[END_REF]:

Ψ w1 (A) = [f(A)] T Υ f(A) (1.34)
A weighed direct-fit criterion can be found in [Yeredor, 2000]. In practice, it could be difficult to obtain an optimal weight matrix Υ since the statistical characterization of the source components is not fully available. Comprehensive discussion on choosing the best weighting matrix for each criterion according to different statistical models of the sources is out of the scope of this thesis report. For more details, one may refer to the following references [Yeredor, 2000, Tichavský and Yeredor, 2009, Slapak and Yeredor, 2011, Yeredor, 2012].

Nonnegative matrix factorization

In the previous section we have addressed the ICA problem. The effectiveness of ICA relies on the statistical independence constraint of the source components. However, in ICA no constraint is made on the sign of the components either of the mixing matrix or of the sources. In other words, the negative values are allowed in the data representation model. Nevertheless, the negative values of both observations and source components are physically meaningless in a wide variety of real-world data. For example, the color and intensity of a pixel in an image is given by a nonnegative number; the frequency of a word that appears in a document is assigned by a nonnegative value; and a music note in a piece of music contributes a nonnegative amount to the signal power spectrum. In addition, some observed signal is most naturally characterized by the exclusively additive combination of the source components, for instance, the eye, nose and mouth of a face image, leading to a nonnegative mixing matrix as well. Imposing the nonnegativity constraint in some BSS problems allows us to obtain results in agreement with physical reality and to facilitate easier interpretation. However, we often lose the nonnegativity constraint during the ICA procedure, for example, when we perform whitening by means of PCA [START_REF] Plumbley | Non-negative mixtures[END_REF].

In this section, we discuss a specific method for the use of nonnegativity constraints in BSS problems, namely Nonnegative Matrix Factorization (NMF). Let us rewrite the linear BSS model (1.1) in the following matrix form:

X = AS + E (1.35)
where X ∈ R N ×M is an observation matrix whose (n, m)-th component denotes the m-th realization of the n-th component of the observation vector x (1.1), where A ∈ R N ×P continues to express the mixing matrix, where S ∈ R P ×M is a source matrix whose 1.4. Nonnegative matrix factorization

Chapter 1

(p, m)-th component gives the value for the corresponding realization of the p-th source, and where the noise matrix E ∈ R N ×M can be similarly defined. We are interested in the conditions where both the mixing matrix A and the source matrix S have nonnegative components. Therefore, we can define the following NMF problem:

Problem 1.3. Given a nonnegative observation matrix X ∈ R N ×M +
, find a mixing matrix A ∈ R N ×P and a source matrix S ∈ R P ×M , such that:

X = AS + E (1.36)
by minimizing the residual term E, subject to A and S having nonnegative components.

NMF was initially investigated by Paatero and Tapper as the concept of positive matrix factorization [Paatero andTapper, 1994, Paatero, 1997], and gained popularity through the well-known works of Lee and Seung [Lee and Seung, 1999]. Various methodologies and extensions of NMF have been constantly emerging during the last few years. The existing NMF algorithms can be divided into the following four categories, including the basic NMF, the constrained NMF, the structured NMF and the generalized NMF [START_REF] Wang | Nonnegative matrix factorization: A comprehensive review[END_REF]. NMF and its various extensions have become prominent tools in a wide range of areas, such as image classification [Das [START_REF] Gupta | [END_REF], face recognition [START_REF] Zafeiriou | Exploiting discriminant information in nonnegative matrix factorization with application to frontal face verification[END_REF], music transcription [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: with application to music analysis[END_REF], spectroscopy separation [ Van Benthem and Keenan, 2004], and pattern recognition [Cichocki et al., 2009]. Many systematic surveys of NMF have been addressed in literature, such as the review papers [Tropp, 2003, Chu et al., 2004, Berry et al., 2007, Buciu, 2008, Plumbley et al., 2010, Wang and Zhang, 2013, Zhou et al., 2014, Smaragdis et al., 2014], and an in-depth book [Cichocki et al., 2009].

In this section, we provide a concise overview of limited areas of NMF which are more relevant to the work described in this thesis, including basic NMF and one generalized NMF, namely semi-NMF. For a more broad discussion one may refer to the aforementioned literature.

Existence and uniqueness

Given a nonnegative matrix X, its trivial NMF solution of form X = AS always exists, for instance A = X and S = I. Vasiloglou et al. proved that every nonnegative matrix has a non-trivial NMF solution by relating NMF to the Completely Positive Matrix (CPM) factorization [START_REF] Vasiloglou | Nonnegative matrix factorization, convexity and isometry[END_REF]. It shows that the set of CPMs forms a convex cone and the solution of NMF belongs to the CPM cone, leading to the conclusion that solving NMF is a convex optimization problem. However, finding a practical description of the CPM cone is still an open question. So although the problem could be reformulated as a convex problem, there is no algorithm known for solving it [START_REF] Vasiloglou | Nonnegative matrix factorization, convexity and isometry[END_REF].

NMF is considered to be essentially unique when its solution is only subject to scale and permutation indeterminacies. Unfortunately, NMF solutions often suffer from additional rotation indeterminacy. For example, there are many ways to select a ratation matrix R which is not necessarily nonnegative or not necessarily equal to a product of a diagonal matrix and a permutation matrix, so that

X = AS = ARR -1 S = A 1 S 1 (1.37)
Chapter 1

1.4. Nonnegative matrix factorization and the rotated matrices A 1 = AR = A and S 1 = R -1 S = S are still nonnegative [Cichocki et al., 2009]. That is to say, only the nonnegativity constraint is not sufficient to guarantee the uniqueness. Nevertheless, incorporating some extra constraints, such as separability, sparsity or symmetry, is sufficient to solve the NMF problem uniquely. The first sufficient condition of uniqueness for separable NMF was provided by [START_REF] Donoho | When does nonnegative matrix factorization give a correct decomposition into parts[END_REF]]. Laurberg et al. presented a stronger necessary and sufficient condition for the uniqueness based on the boundary closeness assumption of factors [START_REF] Laurberg | Theorems on positive data: On the uniqueness of NMF[END_REF]. Gillis analyzed the uniqueness of sparse NMF from a geometric point of view [Gillis, 2012]. More recently Huang et al. derived new uniqueness results for symmetric and asymmetric NMF, which stated that a sufficient condition for uniqueness was that the conic hull of the latent factors was a superset of a particular second-order cone [START_REF] Huang | Nonnegative matrix factorization revisited: Uniqueness and algorithm for symmetric decomposition[END_REF]. In addition checking this condition was shown to be Non-deterministic Polynomial time (NP)-complete. We omit these mathematical details in this report.

Basic NMF algorithms

The idea of basic NMF is to find efficient solutions to the NMF problem under the sole nonnegativity constraint. Due to the NP-hardness and lack of appropriate convex formulation of the NMF problem, the non-convex formulations are adopted in practice.

Generally the non-convex formulations are relatively easy to solve but usually only local minima are achievable. Practical NMF algorithms perform the optimization by minimizing proper cost functions that measure the distance or discrepancy between the given matrix X and the factorization AS. The optimization schemes differ from each other in the cost functions and the optimization procedures.

Cost functions

The Square of Euclidean Distance (SED) in equation (1.38) and the Generalized Kullback-Leibler Divergence (GKLD) in equation (1.39) are the most commonly used cost functions [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]:

Ψ F (A, S) = || X -AS || 2 F = N n=1 M m=1 ( X n,m -[AS] n,m ) 2 (1.38) Ψ KL (A, S) = N n=1 M m=1 ( X n,m log X n,m [AS] n,m -X n,m + [AS] n,m ) (1.39)
where X n,m is the (n, m)-th component of X, and where [AS] n,m denotes the (n, m)-th component of the matrix product AS. Other cost functions include Itakura-Saito Divergence (ISD) [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: with application to music analysis[END_REF], α-divergence [START_REF] Cichocki | Non-negative matrix factorization with α-divergence[END_REF], β-divergence [Kompass, 2007], γ-divergence [START_REF] Cichocki | Families of alpha-betaand gamma-divergences: Flexible and robust measures of similarities[END_REF], Earth mover's distance [START_REF] Sandler | Nonnegative matrix factorization with Earth mover's distance metric for image analysis[END_REF]] and many others. Different cost functions can be chosen according to the prior information of the probability distribution of the noises. For example, the minimizations of the SED, GKLD and ISD are equivalent to the maximum likelihood estimators under Gaussian additive, Poisson and multiplicative Gamma noises, respectively [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF]. Although generally these cost functions are not jointly 1.4. Nonnegative matrix factorization

Chapter 1 convex in both A and S, they are separately convex in either A or S, giving birth to a wide range of NMF algorithms based on an alternating minimization scheme.

Algorithms

Lee and Seung introduced the SED and GKLD cost functions, and applied Multiplicative Updates (MU) iteratively between A and S [Lee and Seung, 1999]. For example, the MU rule for SED is shown as follows:

A (it+1) = A (it) [X(S (it) ) T ] [A (it) S (it) (S (it) ) T ] (1.40) S (it+1) = S (it) [(A (it+1) ) T X] [(A (it+1) ) T A (it+1) S (it) ] (1.41)
where it is the iteration index, where denotes Hadamard product, and where means matrix element-wise division. The MU rules can be seen as adaptive rescaled gradient descent algorithms and they have been widely used as a baseline due to their simplicity.

Their multiplicative natures guarantee that the resulting matrix components cannot become negative. Berry et al. analyzed the convergence of the MU algorithms of Lee and Seung [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF].

Considering the separate convexity of the SED cost function (1.38), the joint optimization can be converted into an alternate minimization of (1.38) with respect to A or B while fixing the other one in the least-square sense under the nonnegativity constraint. Such an alternate optimization scheme is referred to as Alternating Nonnegative Least Squares (ANLS) [START_REF] Kim | Algorithms for nonnegative matrix and tensor factorizations: a unified view based on block coordinate descent framework[END_REF]. The following two-step algorithm is the simplest way to achieve ANLS [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF]:

A (it+1) = P + { X(S (it) ) T [S (it) (S (it) ) T ] -1 } (1.42) S (it+1) = P + { [(A (it+1) ) T A (it+1) ] -1 (A (it+1) ) T X } (1.43)
where P + {•} sets all negative values of the input matrix to zero. Kim et al. claimed that the projection of negative values onto the positive orthant can not guarantee the decrease of the cost function [START_REF] Kim | Fast Newton-type methods for the least squares nonnegative matrix approximation problem[END_REF]. In practice, this procedure was reported to have good performance [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF]. Rather than using the projection operation, the Nonnegativity constrained Least-Square (NLS) problem can be solved by means of the active set algorithm proposed by [START_REF] Lawson | Solving Least Squares Problems[END_REF]. Fast NLS algorithms were derived by Bro and De Jong [Bro and De Jong, 1997] and by Van Benthem and Keenan [Van Benthem and Keenan, 2004], respectively. By adopting the active set based NLS algorithm, Kim and Park proposed a framework of ANLS methods for NMF [Kim and Park, 2008].

In order to improve the performance, as well as to accelerate the convergence rate, many NMF methods were proposed based on sophisticated optimization schemes. Merritt and Zhang proposed an interior-point gradient method for solving the NLS problem [START_REF] Merritt | Interior-point gradient method for large-scale totally nonnegative least squares problems[END_REF]. Gonzalez and Zhang introduced a multiplicative regulatory factor leading to an accelerated MU algorithm [START_REF] Gonzalez | Accelerating the Lee-Seung algorithm for non-negative matrix factorization[END_REF]. Chu et al. derived a family of Newton-type approaches for NMF, based on a projected Newton method, the use of an alternating direction iterative Newton method and sequential quadratic programming, respectively [START_REF] Chu | Optimality computation and interpretation of nonnegative matrix factorizations[END_REF] This differentiable parametrization effectively transforms the nonnegativity constrained optimization problem into a problem with no constraint at all, leading to an alternating gradient descent method [START_REF] Chu | Optimality computation and interpretation of nonnegative matrix factorizations[END_REF]. Lin proposed an ANLS algorithm by means of projected gradient methods with an adaptive gradient decent step-size, which achieved a fast convergence rate [Lin, 2007]. In order to obtain a more accurate solution, Zdunek and Cichocki proposed to use the second-order approximations in the Taylor expansion of the cost function, and proposed a family of projected quasi-Newton optimization approaches [START_REF] Zdunek | Nonnegative matrix factorization with constrained second-order optimization[END_REF]. Li and Zhang developed a fixed-point NMF approach based on an exact block coordinate descent method, namely FastNMF [START_REF] Li | FastNMF: highly efficient monotonic fixed-point nonnegative matrix factorization algorithm with good applicability[END_REF]. Cichocki and Phan proposed a class of Hierarchical Alternating Least Squares (HALS) algorithms, giving birth to FastHALS NMF approaches [Cichocki and Phan, 2009]. It is known that a best rank-one approximation can be achieved in polynomial time. Gillis and Glineur proposed to identify an optimal rank-one NMF solution and subtract it from the input matrix. Then the similar procedure was repeated on the residual matrix under an upper bound constraint [START_REF] Gillis | Using underapproximations for sparse nonnegative matrix factorization[END_REF]. Guan et al. provided a Nesterov's optimal gradient method, in which the step-size was determined by the Lipschitz constant hence avoiding the classical time consuming line search procedure [START_REF] Guan | NeNMF: An optimal gradient method for nonnegative matrix factorization[END_REF]. Moreover in order to corp with large-scale problems, many specific optimization approaches were developed, such as an interior-point method [START_REF] Merritt | Interior-point gradient method for large-scale totally nonnegative least squares problems[END_REF], local optimization [Cichocki and Phan, 2009], parallel computing [START_REF] Liu | Distributed nonnegative matrix factorization for web-scale dyadic data analysis on MapReduce[END_REF], low-rank approximation [START_REF] Zhou | Fast nonnegative matrix/tensor factorization based on low-rank approximation[END_REF], Alternating Direction Method of Multipliers (ADMM) [Xu et al., 2012, Sun andFévotte, 2014], to cite only a few.

Semi-NMF algorithms

Basic NMF restricts every component in the observation matrix X to be nonnegative. When the matrix X is unconstrained, which means it may have mixed signs, Ding et al. suggested to restrict A to be nonnegative while placing no restriction on the signs of S, leading to the following semi-NMF problem [Ding et al., 2010]:

Problem 1.4. Given an observation matrix X ∈ R N ×M , find two matrices A ∈ R N ×P and S ∈ R P ×M , such that

X = AS + E (1.44)
by minimizing the residual term E, subject to A having nonnegative components.

The motivation of semi-NMF was originated from the perspective of clustering with A denoting the cluster indicators and S representing the cluster centers [START_REF] Li | The relationships among various nonnegative matrix factorization methods for clustering[END_REF]. In this case, the SED cost function (1.38) can be viewed as an objective function of K-means clustering. On the other side, semi-NMF is also physically interpretable in some practical applications where the observation data is not always nonnegative, so the source components might also have some negative elements reflecting the phase information [START_REF] Wang | Nonnegative matrix factorization: A comprehensive review[END_REF]. By using the SED cost function (1.38), the semi-NMF algorithm of Ding et al. [Ding et al., 2010] iterates between finding the optimum S for a given A, which is given by the classic least-square solution:

S (it+1) = [(A (it) ) T A (it) ] -1 (A (it) ) T X (1.45)
Chapter 1

and updating the estimate of A for a given S by the MU rule:

A (it+1) = A (it) [X(S (it+1) ) T ] + + A (it) [S (it+1) )(S (it+1) ) T ] - [X(S (it+1) ) T ] -+ A (it) [S (it+1) )(S (it+1) ) T ] + (1.46)
where U + and U -correspond to the positive and negative parts of the input matrix U , respectively, given by:

U + m,n = (|U m,n | + U m,n )/2, U - m,n = (|U m,n | -U m,n )/2 (1.47)
The convergence of this alternating minimization procedure was also proved. Park and Kim treated semi-NMF as a one-sided NMF problem and solved it with the NLS method [START_REF] Park | One-sided non-negative matrix factorization and non-negative centroid dimension reduction for text classification[END_REF]. Wang and Li introduced a random projection technique to accelerate the semi-NMF [START_REF] Wang | Efficient nonnegative matrix factorization with random projections[END_REF]. Kumar et al. proposed a semi-NMF algorithm based on a max-margin classification criterion [START_REF] Kumar | Max-margin semi-NMF[END_REF]. Additionally, semi-NMF finds its applications in template matching [Le Roux et al., 2008], image resolution enhancement [START_REF] Bevilacqua | Neighbor embedding based single-image super-resolution using seminonnegative matrix factorization[END_REF], motion segmentation [START_REF] Mo | Semi-nonnegative matrix factorization for motion segmentation with missing data[END_REF], graph theory [START_REF] Mankad | Discovery of path-important nodes using structured semi-nonnegative matrix factorization[END_REF], and hyperspectral image unmixing [START_REF] Yokoya | Nonlinear unmixing of hyperspectral data using semi-nonnegative matrix factorization[END_REF].

Nonnegative independent component analysis

In order to achieve uniqueness of the semi-NMF solution, some additional constraints are necessary, such as mutual independence, sparsity or semi-orthogonality [START_REF] Plumbley | Non-negative mixtures[END_REF], Cichocki, 2013]. The mutual statistical independence, as well as the nonnegativity assumption of the source components leads to the following Nonnegative ICA (NICA) problem:

Problem 1.5. Given M realizations of an observation vector x ∈ R N , find a mixing matrix A ∈ R N ×P and the M corresponding realizations of a source vector s ∈ R P , such that:

x = As + ε (1.48)
where s has statistically independent and also nonnegative components, and where ε ∈ R N is an additive noise vector independent of s.

Preprocessing and identifiablity

A source component s i is called nonnegative if p(s i < 0) = 0 and will be called wellgrounded, if it has a non-vanishing pdf in any positive neighborhood of zero, such that for any δ > 0 we have p(s i < δ) > 0 [Plumbley, 2002]. We assume that we have a noiseless observation vector x and the source components s i have unit variance. Plumbley proposed to perform a whitening step on x, yielding z = W x, where W ∈ R P ×N is the classical whitening matrix. Then the whitened observation vector z can be written as z = U s, where U ∈ R P ×P is an orthonormal matrix. Therefore, the NICA problem is reduced to finding an orthonormal matrix V ∈ R P ×P , for which s = V z is equal to s up to scale and permutation indeterminacies.

Chapter 1

1.5. Nonnegative independent component analysis Plumbley analyzed the identifiability condition for an NICA problem in [Plumbley, 2002]. Firstly, the nonnegativity of an orthonormal matrix is stated in the following lemma:

Lemma 1.1. Let U (0) ∈ R P ×P be an orthonormal matrix such that U (0) (U (0) ) T = (U (0) ) T U (0) = I P . Then all the elements of U (0) are nonnegative if and only if U (0) is a permutation matrix.

Then based on lemma 1.1, we can obtain the following NICA identifiability condition:

Theorem 1.4. Let s ∈ R P be a random vector of nonnegative and well-grounded independent source components, each with unit variance, and let

s = V z = (V U )s = U (0) s
be an orthonormal rotation of s, where U (0) (U (0) ) T = (U (0) ) T U (0) = I P . Then U (0) is a permutation matrix if and only if s is nonnegative with probability 1.

Many practical algorithms perform optimizations by minimizing a cost function. The following corollary follows:

Corollary 1.1. For the same conditions as theorem 1.4, suppose that Ψ(V ) is a cost function such that Ψ(V ) = 0 if and only if s = V z is nonnegative with probability 1. Then Ψ(V ) = 0 if and only if U (0) = V U is a permutation matrix, such that s is a permutation of the source vector s.

NICA algorithms

In order to derive a suitable cost function, Plumbley suggested to project s onto the positive orthant by s + = [ s + 1 , . . . , s + P ] with s + i = max(0, s i ), and to construct a reestimate of z = V T s by z = V T s + [Plumbley, 2002]. Then the cost function is given as follows:

Ψ NICA (V ) = || z -z || 2 F = z -V T s + 2 F (1.49)
The value of (1.49) will be zero if V is obtained such that all the s i are nonnegative or s = s + . Plumbley proposed to minimize criterion (1.49) by means of axis rotations and geodesic search over the Stiefel manifold of orthogonal matrices [Plumbley, 2003].

The convergence of the geodesic NICA algorithm was later proved in [START_REF] Ye | Monotonic convergence of a nonnegative ICA algorithm on Stiefel manifold[END_REF].

Yuan and Oja extended the FastICA algorithm for solving NICA problem [START_REF] Yuan | A FastICA algorithm for nonnegative independent component analysis[END_REF]. In [START_REF] Plumbley | A nonnegative PCA algorithm for independent component analysis[END_REF], Plumbley and Oja pointed out that the stochastic gradient algorithm for minimizing (1.49) was actually a special case of the nonlinear PCA algorithm which was earlier investigated by Oja [Oja, 1997]. The same authors further showed that in the Stiefel manifold of orthogonal matrices, the cost function had no local minimum and it was a Lyapunov function for the matrix gradien flow. Consequently, they proposed a gradient descent algorithm which monotonically converged to the global minimum [START_REF] Oja | Blind separation of positive sources by globally convergent gradient search[END_REF]. The global convergence of such an algorithm was proved by Ye [Ye, 2006]. In [Plumbley, 2005], Plumbley explored the geometrical methods to tackle the NICA problem, such as the manifold and Lie group of special orthogonal matrices. Winther and Petersen presented an empirical Bayesian framework for NICA 1.6. Semi-nonnegative ICA and nonnegative JDC

Chapter 1 [START_REF] Winther | Bayesian independent component analysis: Variational methods and non-negative decompositions[END_REF]. Ouedraogo et al. developed an axis pair rotation method for NICA by means of Newton optimization [START_REF] Ouedraogo | Non-negative independent component analysis algorithm based on 2D Givens rotations and a Newton optimization[END_REF]. A regularized gradient descend algorithm for NICA was presented in [START_REF] Ouedraogo | Regularized gradient algorithm for non-negative independent component analysis[END_REF]. In many practical applications, the nonnegative source components are not well-grounded, Zheng et al. then derived an NICA algorithm minimizing the mutual information [START_REF] Zheng | Nonnegative independent component analysis based on minimizing mutual information technique[END_REF]. Yu et al. formulated the NICA problem as a maximization problem of kurtosis under the nonnegativity constraint of the estimated sources, and solved it by means of a Lagrange multiplier method [START_REF] Yu | Non-negative independent component analysis and its application[END_REF]. Furthermore, NICA plays important roles in image separation [Plumbley, 2003], music transcription [Plumbley, 2003, Dittmar andUhle, 2004], hyperspectral image demixing [START_REF] Bakir | Non-negative maximum likelihood ICA for blind source separation of images and signals with application to hyperspectral image subpixel demixing[END_REF], gene module identification [START_REF] Gong | Gene module identification from microarray data using nonnegative independent component analysis[END_REF], ElectroEncephalography (EEG) / MagnetoEncephaloGraphy (MEG) source imaging [START_REF] Valdés-Sosa | EEG source imaging with spatio-temporal tomographic nonnegative independent component analysis[END_REF], chemical spectra analysis [START_REF] Shao | Extraction of chemical information from complex analytical signals by a non-negative independent component analysis[END_REF], and remote sensing [START_REF] Yu | Remote sensing image fusion based on integer wavelet transformation and ordered nonnegative independent component analysis[END_REF].

1.6 Semi-nonnegative ICA and nonnegative JDC

The semi-nonnegative ICA problem

In this thesis report, we tackle another important nonnegativity constrained ICA problem. Rather than imposing both nonnegativity and statistical independence on the source vector s as NICA, we propose to exploit the statistical independence of s and the nonnegativity of the mixing matrix A, leading to what we call Semi-Nonnegative ICA (SeNICA) [Coloigner, 2012, Coloigner et al., 2014a]. Similarly to semi-NMF, the term semi refers to that only the mixing matrix A is constrained to have nonnegative components in the model. More precisely, the SeNICA problem is defined as follows:

Problem 1.6. Given M realizations of an observation vector x ∈ R N , find a mixing matrix A ∈ R N ×P and the M corresponding realizations of a source vector s ∈ R P , such that:

x = As + ε (1.50)
where A is constrained to have nonnegative components, where s has statistically independent components, and where ε ∈ R N is an additive noise vector independent of the source vector s.

The differences between SeNICA and ICA, NICA, NMF, semi-NMF mainly consist in the assumptions chosen to solve the BSS problem x = As + ε, which are summarized in table 1.1. The SeNICA problem is often encountered in BSS applications where only the mixing matrix is considered to be nonnegative. For instance, in Magnetic Resonance Spectroscopy (MRS), the column vectors of A represent the positive concentrations of the spectra of the source chemical components, while the source spectra s may have negative amplitudes, for example in the case of J-modulation [Sajda et al., 2004, Horská andTkác, 2012].

Of course, if the SeNICA problem 1.6 is noiseless, and if A and s fully satisfy the conditions listed in section 1.2.1, the identifiability of ICA under these conditions will imply that the explicit nonnegativity constraint on A would be unnecessary. However, in practice, these theoretical conditions can not be completely met due to the model

Statistical independence of s √ √ √
Well-groundedness of s √ Table 1.1: Some different assumptions of NMF, semi-NMF, ICA, NICA and SeNICA for solving the BSS problem x = As + ε.

error, the additional noise and the finite number of available samples. Hence, the complementary nonnegativity information should be explicitly considered in order to obtain results consistent with physical reality [Coloigner, 2012]. Furthermore, SeNICA has already shown its good ability in separating MRS signals [Wang et al., 2013a, Coloigner et al., 2014a] and in classification of 3D dose distribution for prediction of rectal bleeding in prostate cancer radiotherapy [Coloigner et al., 2014b]. We proposed to solve the SeNICA problem by imposing the nonnegativity constraint through cumulant-based ICA algorithms. As discussed in section 1.2, compared with negentropy, likelihood estimation and mutual information, the cumulants avoid the complicated procedure of estimating the pdf and the time consuming integral computation, thus they are easier to compute. Following the study in section 1.3, the properties of cumulants allow us to solve the SeNICA problem by means of nonnegative JDC methods, which will be addressed in the next section.

Nonnegative JDC algorithms

The Nonnegative JDC (NJDC) problem is defined by imposing the nonnegativity constraint on the joint transformation matrix, as follows:

Problem 1.7. Given K symmetric matrices C (k) ∈ R N ×N and an integer P , find a joint transformation matrix A ∈ R N ×P and K diagonal matrices D (k) ∈ R P ×P such that for each index k ∈ {1, 2, • • • , K}, we have:

C (k) = AD (k) A T (1.51)
subject to A having nonnegative components.

The NJDC problem can be solved by minimizing the criterion functions presented in section 1.3 with additional nonnegativity constraint. The existing NJDC algorithms are mainly contributed by the french research team of Laboratoire Traitement du Signal et de l'Image (LTSI) in Rennes. Coloigner et al. resorted to use the following two schemes in order to impose nonnegativity on A in the direct-fit JDC criterion (1.29) [Coloigner, 2012, Coloigner et al., 2014a, Coloigner et al., 2014c]: variable has been originally adopted by Chu et al. for solving NMF problem [START_REF] Chu | Optimality computation and interpretation of nonnegative matrix factorizations[END_REF]. In addition, its effectiveness in the context of nonnegative CPD (without symmetry constraints) was shown by Royer et al. [START_REF] Royer | Computing the polyadic decomposition of nonnegative third order tensors[END_REF], Royer, 2013]. Coloigner et al. proposed to minimize the direct-fit criterion (1.29) with square and exponential changes of variable, respectively, by applying two fundamental optimization strategies: line search and trust region procedures. Regarding the linear search method, they developed gradient descent, conjugate gradient and Newton-like algorithms combined with a search for optimal step-sizes [Coloigner et al., 2014c]. Moreover, two Alternating Least Squares (ALS) procedures also accelerated with an Enhanced Line Search (ELS) were derived [START_REF] Coloigner | Semi-nonnegative independent component analysis: the (3,4)-SENICA exp method[END_REF][START_REF] Coloigner | Décomposition de tableaux d'ordre trois semi-nonnégatifs et semi-symétriques[END_REF], Coloigner et al., 2011a, Coloigner et al., 2014a]. Regarding the trust region approach, Coloigner et al. proposed a Levenberg Marquardt (LM) method based on an approximation of the Hessian inside a trust region [START_REF] Coloigner | Décomposition de tableaux d'ordre trois semi-nonnégatifs et semi-symétriques[END_REF],Coloigner et al., 2014c]. Nevertheless, like all the methods based on line search and trust region strategies, those methods appear to be dependent on initialization, and they therefore require a multi-initialization procedure in practice, leading to an increase of numerical complexity [Coloigner, 2012].

A = B B = B 2 , A = exp(B) (1.

Contribution and outline of thesis

In this thesis report, we also adopt the square change of variable A = B 2 due to its simplicity. Then the NJDC problem 1.7 can be transformed into the following unconstrained one:

Problem 1.8. Given K symmetric matrices C (k) ∈ R N ×N and an integer P , find a matrix B ∈ R N ×P and K diagonal matrices D (k) ∈ R P ×P such that for each index k ∈ {1, 2, • • • , K}, we have:

C (k) = B 2 D (k) (B 2 ) T (1.53)
Then the joint transformation matrix A = B 2 is nonnegative for sure.

In order to overcome the drawback of existing NJDC methods, we resort to use Jacobilike optimization procedures, which are known to be less sensitive to initialization. More precisely, the matrix B is decomposed as a product of elementary matrices:

B = i,j Θ (i,j) (ξ i,j ) (1.54)
where Θ (i,j) (ξ i,j ) is a unit-determinant elementary matrix determined by only one parameter ξ i,j , and where (i, j) is the pivot index. The decomposition (1.54) can be achieved by polar, LU or QR matrix factorizations. In a Jacobi-like optimization procedure, each elementary matrix Θ (i,j) (ξ i,j ) is estimated sequentially by minimizing a proper criterion function. In fact, in our case, the parameter of one elementary matrix can be simply estimated by searching the roots of a polynomial function. We proposed five Jacobi-like methods based on different factorizations of B incorporating different cost functions defined in section 1.3, namely JDC + LU-3 [Wang et al., 2013a], iJDC + LU-1 [Wang et al., 2014a], iJDC + QR-1 [START_REF] Wang | CP decomposition of semi-nonnegative semi-symmetric tensors based on QR matrix factorization[END_REF], Wang et al., 2014a], JDC + LU-1 [START_REF] Wang | Decomposition of semi-nonnegative semi-symmetric three-way tensors based on LU matrix factorization[END_REF][START_REF] Wang | Parameterization of the inverse of the nonnegative matrix for nonnegative joint diagonalization[END_REF], and JDC + LU-2 [START_REF] Wang | Parameterization of the inverse of the nonnegative matrix for nonnegative joint diagonalization[END_REF]. Besides those five Jacobi-like methods

Chapter 1

1.6. Semi-nonnegative ICA and nonnegative JDC based on a square change of variable, we developed a sixth algorithm, based on the Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], that tackles the NJDC problem 1.7 by minimizing the direct-fit JDC criterion (1.29). Such a method, namely JDC + ADMM-3 , estimates A directly, and is well suited for large-scale problems. The differences between these six algorithms are summarized in table 1.2. The names of the algorithms are straightforwardly determined by the cost functions, optimization schemes and preprocessing step that they used. Taking the algorithm <iJDC + LU-1 > for example, i) the prefix "i" means that it requires to inverse each input matrix C (k) as a preprocessing step; ii) the superscript "+" represents the nonnegativity constraint; iii) the subscript "LU" denotes that the square change of variable A = B 2 is used and B is decomposed using the LU matrix factorization; iv) the ensuing subscript "1" indicates the first indirect-fit criterion (1.27), while "2" and "3" signify the second indirect-fit criterion (1.28) and the direct-fit criterion (1.29), respectively. Furthermore, in some practical SeNICA problems, the dimension of the observation space must be reduced. However, the classical dimension compression procedure, such as whitening by PCA, breaks the nonnegativity property of the compressed mixing matrix. We proposed a new nonnegative compression method in [START_REF] Wang | Nonnegative compression for semi-nonnegative independent component analysis[END_REF], which guarantees the nonnegativity of the compressed mixing matrix. This thesis report is organized as follows:

• Chapter 2: This chapter is devoted to the detailed description of the proposed NJDC methods. First of all, we introduce some necessary preliminaries for our methods, including LU and QR matrix factorizations, the Jacobi-like optimization principle and the concept of ADMM. Next, the proposed algorithms are presented. We separate the six algorithms into four groups for the purpose of convenience.

For each group, we provide the reformulated cost functions and the comprehensive derivation of the algorithms. The assumptions on each group of algorithms as well as the practical issues are also addressed. A numerical complexity study of the methods proposed in this thesis is given at the end of this chapter.

• Chapter 3: Our algorithms are evaluated on random synthetic matrices. Several classical nonorthogonal JDC methods without nonnegativity constraint, and one NJDC method proposed by Coloigner et al. [Coloigner et al., 2014c], are tested as baseline algorithms. The performance is assessed in terms of the matrix estimation accuracy and the numerical complexity. The convergence property, the influence of Signal to Noise Ratio (SNR), the impact of dimensions, the effect of coherence, and the influence of the modulus of uniqueness, are extensively studied through Monte Carlo experiments. The obtained results show that the proposed algorithms offer better estimation accuracy by exploiting the nonnegativity a priori, especially in difficult contexts.

• Chapter 4: The efficiency of SeNICA methods is illustrated by solving three BSS problems from very different application domains. The first one which is rooted in the biomedical diagnostics, consists of separating synthetic MRS spectra. The second one is an automatic music transcription problem. It aims at estimating the musical notes and their attributes from a single-channel music record by decomposing its magnitude spectrogram. The third one is devoted to separating No is required.

iJDC + LU-1 √ √ √ iJDC + QR-1 √ √ √ JDC + LU-1 √ √ √ JDC + LU-2 √ √ √ JDC + LU-3 √ √ √ JDC + ADMM-3 √ √ √
Chapter 1 1.7. List of publications digital images which are degraded by the so-called show-through effect. Such a problem is usually caused by the seeping of ink from the reverse side or scanning a semi-transparent paper. The behavior of the proposed methods are compared with classical ICA, NICA, NMF and semi-NMF methods. The experimental results show that by fully exploiting the prior information of data, such as the nonnegativity of the mixing matrix and the statistical independence of the sources, the proposed methods achieve better estimation results.

• Chapter 5: This final section concludes the thesis report by providing a closing discussion on the contributions presented in the chapters as well as possible future extensions.

1.7 List of publications Chapter 2

Nonnegative joint diagonalization by congruence

This chapter is devoted to provide an in-depth scrutiny of the proposed six novel NJDC algorithms. The first five algorithms impose the nonnegativity constraint of the joint transformation matrix A by means of a square change of variable A = B 2 and estimate B using a Jacobi-like optimization procedure. They differ in the cost functions, matrix factorization schemes and preprocessing steps. The sixth algorithm estimates A directly through ADMM, which is an improved variant of the classical augmented Lagrangian method.

The remainder of this chapter is organized as follows. In section 2.1, we introduce some necessary preliminaries for the proposed methods, including LU and QR matrix factorizations, the Jacobi-like optimization principle and the concept of ADMM. Then the six algorithms are presented in the following order, JDC + LU-3 in section 2.2, iJDC +

LU-1

and iJDC + QR-1 in section 2.3, JDC + LU-1 and JDC + LU-2 in section 2.4, and JDC + ADMM-3 in section 2.5. In each section, firstly we present the comprehensive derivation of the algorithms for estimating a square joint transformation matrix A ∈ R N ×N + . Secondly we discuss the possible generalization of the algorithms to the non-square case A ∈ R N ×P + with N > P , as well as some practical issues. In addition, we provide a numerical complexity analysis of the methods proposed in this thesis in section 2.6.

Preliminaries

LU and QR matrix factorizations

Let us recall the following elementary unit-determinant matrices:

Definition 2.1. A unit upper (or lower) triangular matrix is a upper (or lower, respectively) triangular matrix whose main diagonal elements are equal to 1.

Definition 2.2. An elementary upper triangular matrix with parameters {i, j, u i,j } and i < j is a unit upper triangular matrix whose non-diagonal elements are zeros except the (i, j)-th entry, which is equal to u i,j .

Chapter 2 U (i,j) (u i,j ) with i < j and i, j ∈ {1, 2, . . . , N } denotes an elementary upper triangular matrix:

U (i,j) (u i,j ) =          I i-1 . . . 0 . . . 0 . . . 1 . . . u i,j . . . 0 . . . I j-i-1 . . . 0 . . . 0 . . . 1 . . . 0 . . . 0 . . . I N -j          (2.1)
Definition 2.3. An elementary lower triangular matrix with parameters {i, j, i,j } and i > j is a unit lower triangular matrix whose non-diagonal elements are zeros except the (i, j)-th entry, which is equal to i,j .

L (i,j) ( i,j ) with j < i and i, j ∈ {1, 2, . . . , N } corresponds to an elementary lower triangular matrix:

L (i,j) ( i,j ) =          I i-1 . . . 0 . . . 0 . . . 1 . . . 0 . . . 0 . . . I j-i-1 . . . 0 . . . i,j . . . 1 . . . 0 . . . 0 . . . I N -j          (2.2)
Definition 2.4. A Givens rotation matrix with parameters {i, j, θ i,j } and i < j is equal to an identity matrix except for the (i, i)-th, (j, j)-th, (i, j)-th and (j, i)-th entries, which are equal to cos(θ i,j ), cos(θ i,j ), -sin(θ i,j ) and sin(θ i,j ), respectively.

Q (i,j) (θ i,j ) with i < j and i, j ∈ {1, 2, . . . , N } indicates the corresponding Givens rotation matrix:

Q (i,j) (θ i,j ) =          I i-1 . . . 0 . . . 0 . . . cos(θ i,j )
. . .

-sin(θ i,j ) . . .

0 . . . I j-i-1 . . . 0 . . . sin(θ i,j ) . . . cos(θ i,j ) . . . 0 . . . 0 . . . I N -j          (2.3)
Now, let us consider the following lemmas:

Lemma 2.1. Any (N × N ) unit upper (or lower) triangular matrix can be factorized as a product of N (N -1)/2 elementary upper (or lower, respectively) triangular matrices.

The proof of lemma 2.1 is straightforward by reducing the unit upper (or lower) triangular matrix into an identity matrix by means of Gaussian elimination [Golub and Van Loan, 1996, Chapter 3].

Lemma 2.2. Any (N × N ) orthonormal matrix can be factorized as a product of, at most, N (N -1)/2 Givens rotation matrices.

The proof of lemma 2.2 can be found in [Vaidyanathan, 1993, Chapter 14]. Furthermore, the above two decompositions are not unique.

Chapter 2

Preliminaries

Then we review the following two classical matrix factorization methods for a nonsingular matrix B ∈ R N ×N [Golub and Van Loan, 1996],

LU matrix factorization:

B = LU Λ 1 Π 1 (2.4)
where L ∈ R N ×N is a unit lower triangular matrix, where U ∈ R N ×N is a unit upper triangular matrix, where Λ 1 ∈ R N ×N is a diagonal matrix, and where Π 1 ∈ R N ×N is a permutation matrix.

QR matrix factorization:

B = QR Λ 2 Π 2 (2.5)
where Q ∈ R N ×N is an orthonormal matrix, where R ∈ R N ×N is a unit upper triangular matrix, where Λ 2 ∈ R N ×N is a diagonal matrix, and where Π 2 ∈ R N ×N is a permutation matrix.

The LU and QR matrix decompositions are important linear algebra tools that are widely used in both scientific and engineering applications.

C (k) = (LU ) 2 (Λ 2 1 Π 1 D (k) Π T 1 Λ 2 1 )[(LU ) 2 ] T = (LU ) 2 D (k) 1 [(LU ) 2 ] T (2.6) C (k) = (QR) 2 (Λ 2 2 Π 2 D (k) Π T 2 Λ 2 2 )[(QR) 2 ] T = (QR) 2 D (k) 2 [(QR) 2 ] T (2.7) where: D (k) 1 = Λ 2 1 Π 1 D (k) Π T 1 Λ 2 1 , D (k) 2 = Λ 2 2 Π 2 D (k) Π T 2 Λ 2 2 (2.8)
It is straightforward to verify that D 2 are also diagonal matrices. Therefore, due to the scale and permutation indeterminacies of the JDC problem, the Hadamard square root of the joint transformation matrix, say B, can be expressed as B = LU and B = QR without loss of generality. Moreover, by incorporating lemma 2.1 and lemma 2.2, we obtain the following two elementary factorizations of B:

B = N j=1 N i=j+1 L (i,j) ( i,j ) N i=1 N j=i+1 U (i,j) (u i,j )
(2.9)

B = N i=1 N j=i+1 Q (i,j) (θ i,j ) N i=1 N j=i+1 U (i,j) (u i,j ) (2.10)
where

L (i,j) ( i,j ) ∈ R N ×N , U (i,j) (u i,j ) ∈ R N ×N and Q (i,j) (θ i,j
) ∈ R N ×N are an elementary lower triangular matrix, an elementary upper triangular matrix and a Givens rotation matrix, respectively. As a consequence, estimating B for the NJDC problem 1.8 is converted to the estimate of N (N -1) parameters: i,j and u i,j for the LU decomposition (2.9), or θ i,j and u i,j for the QR decomposition (2.10). It is noteworthy that by using such elementary decompositions, the trivial zero solution of the NJDC problem is avoided, since B is implicitly constrained to have a unit determinant. Instead of simultaneously computing these N (N -1) parameters, we propose to adopt a Jacobi-like procedure which performs N (N -1) sequential optimizations.

Preliminaries

Chapter 2

A Jacobi-like optimization procedure

Most NJDC algorithms presented in this chapter are Jacobi-like. Therefore we highlight in this section the principle of a Jacobi-like optimization procedure.

The purpose of the NJDC problem 1.8 is to find the Hadamard square root B of the nonnegative joint transformation matrix A. We aim at estimating the optimal matrix B by minimizing a proper criterion function Ψ with respect to B. Ψ(B) could be reformulated from the JDC criteria presented in section 1.3.3 incorporating the square change of variable A = B 2 . According to the discussion in section 2.1.1, B can be decomposed as a product of elementary matrices:

B = i,j Θ (i,j) (ξ i,j ) (2.11)
where Θ (i,j) (ξ i,j ) can be elementary triangular matrices or Given rotation matrices. Instead of computing these elementary matrices simultaneously, a Jacobi-like optimization procedure consists of estimating these elementary matrices successively and building up B gradually as the product of all the estimated elementary matrices [ De Lathauwer, 2010]. More precisely, we can start by an initialization of B, namely B (0) , where the superscript (0) denotes the iteration number. At the it-th Jacobi-like iteration, we estimate one of the elementary matrices, namely Θ (i,j) (ξ i,j ), with a specific (i, j) index couple. The optimal parameter ξ i,j is chosen as the global minimum of a local cost function as follows:

ξ opt i,j = argmin ξ i,j Ψ(B (it-1) Θ (i,j) (ξ i,j )) (2.12)
where B (it-1) is the estimate of B at the previous iteration thus it is a constant term.

Compared to the minimization of the global criterion Ψ(B) with respect to B directly, generally the local optimization (2.12) is much easier to solve. Once ξ opt i,j is obtained, B (it-1) is updated as follows:

B (it) = B (it-1) Θ (i,j) (ξ opt i,j ) (2.13)
Then the (it + 1)-th Jacobi-like iteration continues to compute Θ (i,j) (ξ i,j ) with the next (i, j) index. Such a procedure is repeated until the global criterion Ψ(B) is minimized.

In the classical Jacobi-like iteration, the choice of the order of the index (i, j) is crucial. From the standpoint that the decrease of the cost function shall be as large as possible in each Jacobi-like iteration, it makes sense to choose the pair (i , j ) for which

Ψ(B (it-1) Θ (i ,j ) (ξ i ,j )
) is minimal among all potential (i, j) indices. However, finding the optimal (i, j) index in each iteration is computationally expensive. Therefore, in practice, when going through the different subproblems (2.12), one rather follows a fixed (i, j) sequence, for example:

(1, 2), (1, 3), . . . , (1, N ), (2, 1), (2, 3), . . . , (2, N ), . . . ,

(N -1, N ) (2.

14)

This technique is called the cyclic Jacobi-like algorithm [De Lathauwer, 2010]. The pseudo-code of the cyclic Jacobi-like optimization procedure for a NJDC problem is outlined in the appendix C.1.

Chapter 2

2.1. Preliminaries

The alternating direction method of multipliers

In a finite-dimensional setting, ADMM [Glowinski and Marroco, 1975,Gabay and Mercier, 1976,Boyd et al., 2010] is a method that proposes to combine the decomposability of the dual ascent method [Shor, 1985, Bertsekas, 1999, Nedić and Ozdaglar, 2010] with the superior convergence properties of the method of multipliers [Hestenes, 1969, Miele et al., 1972,Bertsekas, 1996]. The algorithm solves structured convex problems of the following form: minimize

x,y f(x) + g(y)

subject to P x + Qy = c (2.15)
where f(•) and g(•) are two convex functions defined on closed convex subsets of real numbers of dimensions N 1 and N 2 , respectively, with

x ∈ R N 1 , y ∈ R N 2 , P ∈ R N 3 ×N 1 , Q ∈ R N 3 ×N 2 and c ∈ R N 3 .
The optimal value of the problem (2.15) will be denoted as follows:

p * = inf{f(x) + g(y)|P x + Qy = c} (2.16)
where inf{•} denotes infimum. The augmented Lagrangian function of (2.15) is formed as follows:

L ρ (x, y, λ) = f(x) + g(y) + λ T (P x + Qy -c) + ρ 2 ||P x + Qy -c|| 2 F (2.17)
where λ is a Lagrangian multiplier vector and ρ > 0 is a penalty parameter. ADMM consists of the iterations:

x (it) = argmin x L ρ ( x, y (it-1) , λ (it-1) ) (2.18) y (it) = argmin y L ρ ( x (it) , y, λ (it-1) ) (2.19) λ (it) = λ (it-1) + ρ( P x (it) + Qy (it) -c ) (2.20)
The ADMM algorithm resembles the dual ascent method and the method of multipliers [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. It includes an x-minimization step (2.18), a y-minimization step (2.19), and a dual variable update step (2.20). It is noteworthy that the classical Augmented Lagrangian Multiplier (ALM) method [Hestenes, 1969, Rockafellar, 1973, Fortin and Glowinski, 1983] minimizes (2.17) jointly with respect to both x and y. For example, ALM replaces the steps (2.18) and (2.19) by

[x (it) , y (it) ] = argmin x,y L ρ ( x, y, λ (it-1) ) (2.21)
which could be more difficult to solve since it involves both f(x) and g(y). On the other side, in ADMM, x and y are updated in an alternating manner, which specifies the term "alternating direction". ADMM can be regarded as a special version of the method of multipliers where a single Gauss-Seidel pass over x and y is used instead of the usual joint minimization [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. The convergence analysis for ADMM was addressed in [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] and references therein.

2.2. The JDC + LU-3 algorithm Chapter 2

In the following three sections, we present five NJDC algorithms for solving problem 1.8 by means of a square change of variable and a Jacobi-like optimization procedure. For the sake of convenience, firstly we assume that the joint transformation matrix A ∈ R N ×P + in problem 1.8 is a square matrix such that N = P . Given K symmetric matrices

C (k) ∈ R N ×N ,
each algorithm is derived for the purpose of estimating A ∈ R N ×N + . Secondly, we discuss the possible extensions of these algorithms in order to seek a nonsquare matrix A ∈ R N ×P + where N > P .

The JDC +

LU-3 algorithm

Algorithm derivation

Let us impose the nonnegativity constraint in the JDC criteria presented in section 1.3.3. The indirect-fit criteria (1.27) and (1.28) seem difficult to handle since it is quite complicated to impose the nonnegativity constraint on A, meanwhile to perform the optimization on its inverse A -1 . In other words, it is difficult to parameterize the set of matrices which are the inverse of nonnegative matrices. On the other side, the direct-fit JDC criterion (1.29) carries out the optimization on A directly. Therein the nonnegativity of A can be enforced straightforwardly. Therefore, the first and natural idea is to incorporate the square change of variable in the direct-fit JDC criterion (1.29), leading to the following nonnegativity constrained direct-ft JDC criterion:

Ψ + 3 (B, {D (k) }) = K k=1 C (k) -(B 2 )D (k) (B 2 ) T 2 F (2.22)
Minimizing (2.22) with respect to both B and the set {D (k) } is difficult. Consequently, as Yeredor's ACDC algorithm [Yeredor, 2002], the proposed JDC + LU-3 algorithm performs the minimization of (2.22) alternately with respect to B and {D (k) }. In order to estimate B, by means of the LU matrix factorization (2.9) and the Jacobi-like optimization procedure, the high dimensional optimization is reduced to search a sequence of sparse triangular matrices. Therefore, JDC + LU-3 consists of the following two minimization schemes:

1. The AC (Alternating Columns) phase minimizes (2.22) with respect to B via Jacobi-like iterations. In each iteration, one elementary triangular matrix factor L (i,j) ( i,j ) or U (i,j) (u i,j ) of B is identified, while keeping its other factors as well as the diagonal matrices {D (k) } fixed. For example,

(it) i,j = argmin i,j Ψ + 3 ( B (it-1) L (i,j) ( i,j ), {D (k,it-1) } ) (2.23)
The sequential updates of B by right multiplying it by L (i,j) ( i,j ) or U (i,j) (u i,j ) with a different index couple (i, j) affects the j-th column vector of B sequentially. Thus we retain the terminology "Alternating Columns phase" from the ACDC algorithm.

2. The DC (Diagonal Centers) phase minimizes (2.22) with respect to the diagonal matrix set {D (k) } while keeping B fixed:

{D (k,it) } = argmin {D (k) } Ψ + 3 (B (it-1) , {D (k) }) (2.24)
Chapter 2

2.2. The JDC + LU-3 algorithm Hence, JDC + LU-3 is an alternating optimization approach including a Jacobi-like optimization phase.

AC phase

In this phase, we minimize (2.22) with respect to L (i,j) ( i,j ) and U (i,j) (u i,j ). In order to simplify the notations, we define an elementary triangular matrix T (i,j) (t i,j ) such that:

T (i,j) (t i,j ) = L (i,j) ( i,j ), if i > j U (i,j) (u i,j ), if i < j (2.25)
Consequently, the LU decomposition (2.9) of B can be simply denoted as follows:

B = N (N -1) i =j T (i,j) (t i,j ) (2.26)
Instead of simultaneously computing these N (N -1) matrices, each elementary triangular matrix with a selected (i, j) index is identified sequentially. For the sake of convenience, we use A, B and D (k) to denote the estimate of A, B and D (k) at (it -1)-th iteration instead of A (it-1) , B (it-1) and D (k,it-1) , respectively. At the current iteration, in other words, at the it-th iteration, the update of B by one T (i,j) (t i,j ), denoted by B (new) , is defined as follows:

B (new) = BT (i,j) (t i,j ) (2.27)
Inserting (2.27) into the optimization (2.23), the estimation of t i,j only consists of minimizing the following local cost function:

Ψ + 3,AC (t i,j ) = K k=1 C (k) -[( B (new) ) 2 ] D (k) [( B (new) ) 2 ] T 2 F (2.28)
The expression of the Hadamard square of the matrix B (new) is shown in the following proposition.

Proposition 2.1. A (new) = ( B (new) ) 2 = ( BT (i,j) (t i,j )) 2 can be expressed as a function of t i,j as follows:

A (new) = ( B (new) ) 2 = B 2 T (i,j) (t 2 i,j ) + 2 t i,j ( b i b j )e T j (2.29)
where b i and b j denote the i-th and j-th column vectors of B, respectively, and where e j is the j-th column vector of the identity matrix I N .

The proof of proposition 2.1 is provided in appendix A.1. Therefore, the cost function (2.28) can be rewritten as an explicit function of t i,j as follows:

Ψ + 3,AC (t i,j ) = K k=1 F (k,4) t 4 i,j + F (k,3) t 3 i,j + F (k,2) t 2 i,j + F (k,1) t i,j + F (k,0) 2 F (2.30)
2.2. The JDC + LU-3 algorithm Chapter 2

where:

F (k,4) = -d (k) j ( b 2 i )( b 2 i ) T (2.31) F (k,3) = -2 d (k) j ( b 2 i )( b i b j ) T + ( b i b j )( b 2 i ) T
(2.32)

F (k,2) = -d ( k) j ( b 2 i )( b 2 j ) T + ( b 2 j )( b 2 i ) T + 4( b i b j )( b i b j ) T (2.33) F (k,1) = -2 d (k) j ( b 2 j )( b i b j ) T + ( b i b j )( b 2 j ) T (2.34) F (k,0) = C (k) -( B 2 )D (k) ( B 2 ) T (2.35)
and where d

(k) j
denotes the (j, j)-th entry of D (k) . (2.30) can be expressed in the following compact matrix form:

Ψ + 3,AC (t i,j ) = K k=1 F (k) τ i,j 2 F = τ T i,j G τ i,j (2.36)
where:

G = K k=1 (F (k) ) T F (k) (2.37)
is a (5 × 5) symmetric coefficient matrix, and where F (k) is a (N 2 × 5) matrix defined as follows:

F (k) = vec(F (k,4) ), vec(F (k,3) ), vec(F (k,2) ), vec(F (k,1) ), vec(F (k,0) ) (2.38)
in which vec(.) reshapes a matrix into a column vector by stacking its columns vertically. τ i,j is a 5-dimensional parameter vector defined as follows:

τ i,j = t 4 i,j , t 3 i,j , t 2 i,j , t i,j , 1 T (2.39)
Cost function (2.36) shows that Ψ + 3,AC (t i,j ) is a 8-th degree polynomial in t i,j . The global minimum t i,j can be obtained by computing the roots of its derivative and selecting the root yielding the smallest value of (2.36). Once the optimal t i,j is computed, the matrix B (new) is computed using equation (2.27), and A (new) is given by A (new) = ( B (new) ) 2 . Then the AC phase is repeated to estimate T (i,j) (t i,j ) with the next (i, j) index.

DC phase

In this phase, the minimization of (2.24) with respect to the diagonal matrices {D (k) } can be separated into K distinct linear least square subproblems:

Ψ (k) 3,DC (D (k) ) = C (k) -AD (k) A T 2 F , k ∈ {1, 2, . . . , K} (2.40)
The optimal solution of D (k,new) given by Yeredor in [Yeredor, 2002] is:

D (k,new) = Diag ( A A) c (k) (2.41)
Chapter 2

2.2. The JDC + LU-3 algorithm where Diag(•) returns a diagonal matrix whose diagonal elements are given by the input vector, where c (k) = vec(C (k) ), and where A is the current estimate of the joint transformation matrix A.

In this context, one full iteration is referred to the combination of one full AC phase and one DC phase. In practice, several iterations are necessary to ensure convergence. One can stop the algorithm when the value of the global cost function (2.22), or the decrease of (2.22) between two successive iterations falls below a fixed positive threshold. Such a stopping criterion is guaranteed to be met since (2.22) is non-increasing in each sweep. We observed empirically that the proposed algorithm converges linearly.

Generalization to the non-square case

In a JDC problem, a non-square joint transformation matrix A ∈ R N ×P + and K diagonal matrices D (k) ∈ R P ×P with N > P are frequently encountered. Let us recall that A = B 2 . Based on the LU decomposition of B ∈ R N ×P , which yields B = LU , where L ∈ R N ×P is a unit lower triangular matrix, and where U ∈ R P ×P is a unit upper triangular matrix, the JDC + LU-3 algorithm can handle non-square matrices directly. Theoretically, L and U can be decomposed as follows:

L = P j=1 N i=j+1 L (i,j) ( i,j ), U = P i=1 P j=i+1 U (i,j) (u i,j ) (2.42)
where the sizes of all the L (i,j) ( i,j ) and U (i,j) (u i,j ) matrices are (N × N ) and (P × P ), respectively, except that of L (N,P ) ( N,P ) which is (N × P ). Below we give an example:   1.0000 0.1176 0.3750 1.0441 0.5000 1.8966

  =   1 0 0 0.3750 1 0 0 0 1     1 0 0 0 1 0 0.5 0 1     1 0 0 1 0 1.8378   1 0.1176 0 1
Indeed, we can treat all the L (i,j) ( i,j ) and U (i,j) (u i,j ) matrices as (N × N ) matrices, and choose the first P column vectors of their product as the final result. In the above example, we actually estimate four (3 × 3) elementary triangular matrices and choose the first two column vectors of the product as the result:   1.0000 0.1176 0 0.3750 1.0441 0 0.5000 1.8966 1

  =   1 0 0 0.3750 1 0 0 0 1     1 0 0 0 1 0 0.5 0 1     1 0 0 0 1 0 0 1.8378 1     1 0.1176 0 0 1 0 0 0 1  

Practical issues

In practice, we observe that if each input matrix C (k) is almost exactly jointly diagonalizable due to a high SNR, the classical non-constrained JDC methods can also give a nonnegative A with high probability. In this situation, imposing the nonnegativity constraint explicitly is unnecessary and increases the computational burden. Therefore, we propose to relax the nonnegativity constraint by directly decomposing A into elementary LU form instead of using the decompositions of B as follows:

A = N j=1 N i=j+1 L (i,j) ( i,j ) N i=1 N j=i+1 U (i,j) (u i,j ) def = N (N -1) i =j T (i,j) (t i,j ) (2.43)
2.3. The iJDC + LU-1 and iJDC + QR-1 algorithms Chapter 2

By inserting (2.43) into the original direct-fit JDC criterion (1.29), each parameter t i,j can be found by minimizing the following 4-th degree polynomial:

Ψ 3,AC (t i,j ) = K k=1 F (k,4) t 2 i,j + (F (k,2) + 4 d (k) j ( b i b j )( b i b j ) T )t i,j + F (k,0) 2 F (2.44)
The derivation of the compact matrix form of (2.44) is omitted in this report. In practice, it is suggested to compute t i,j by minimizing (2.44) first. If all the elements in the j-th column vector of AT (i,j) (t i,j ) have the same sign ς, the j-th column vector of A (new) is replaced by ς( a j +t i,j a i ), while the other column vectors remain unchanged, where a i and a j are the i-th and j-th column vectors of A, respectively. Otherwise, t i,j is computed by minimizing (2.36). A is then updated using proposition 2.1. The pseudo-code of JDC + LU-3 is given in the appendix C.2 .

The iJDC +

LU-1 and iJDC + QR-1 algorithms

Algorithm derivation

The JDC + LU-3 algorithm is conceptually simple. However, the alternating nature of such a method leads to some drawbacks. For example, its performance appears to depend on initialization, therefore in practice it requires a multi-initialization procedure, leading to an increase of numerical complexity. Moreover, similar to the ACDC algorithm, JDC + LU-3 converges linearly, and occasionally its convergence can be slow. In fact, in many practical ICA problems, the estimate of these diagonal matrices D (k) is not necessary. Therefore, without estimating D (k) , we propose to use the indirect-fit criterion (1.27) in order to derive a pure Jacobi-like NJDC algorithm, inspired by the fast convergence property of this kind of procedure. Furthermore, we would like to use both LU and QR matrix factorizations for the Jacobi-like procedure in order to study the influence of the used matrix factorization. We also would like to avoid the inverse operator in criterion (1.27). Now let us consider the structure of C (k) = AD (k) A T with the following two assumptions:

1. A ∈ R N ×N + is nonsingular; 2. Each D (k) ∈ R N ×N
is nonsingular which means that its main diagonal does not contain any zero entry.

Then each matrix C (k) is nonsingular and its inverse can be expressed as follows:

C (k) -1 = A -T D (k) -1 A -1 (2.45)
We use C (k,-1) to denote (C (k) ) -1 for simplicity. Equation (2.45) shows that C (k,-1) also preserves the jointly diagonalizable structure. Furthermore, instead of A -1 , A serves as the joint diagonalizer. As a consequence, by means of a preprocessing step that inverses each input matrices C (k) , A can be estimated by minimizing the following criterion which is modified from (1.27):

Ψ 1 (A) = K k=1 off(A T C (k,-1) A) 2 F (2.
Ψ + 1 (B) = K k=1 off (B 2 ) T C (k,-1) B 2 2 F (2.47)
By means of the LU and QR matrix factorizations of B, the minimization of (2.47) with respect to B is converted to the estimate of the following two sets of N (N -1) parameters, respectively, 1. i,j and u i,j for the LU decomposition (2.9); 2. θ i,j and u i,j for the QR decomposition (2.10).

Instead of simultaneously computing these parameters, we propose two Jacobi-like procedures which perform a series of N (N -1) sequential optimizations. This yields two new algorithms: i) the first algorithm, named iJDC + LU-1 , estimates each i,j and u i,j successively, and ii) the second one, called iJDC + QR-1 , estimates each θ i,j and u i,j sequentially. Now, the difficulty is how to estimate these of parameters, namely L (i,j) ( i,j ) and U (i,j) (u i,j ) for iJDC + LU-1 , and Q (i,j) (θ i,j ) and U (i,j) (u i,j ) for iJDC + QR-1 . Two points should be noted here:

1. For the iJDC + LU-1 algorithm, L (i,j) ( i,j ) and U (i,j) (u i,j ) belong to the same category of matrices, therefore they can be estimated by the same algorithmic procedure just with an emphasis on the relation between the i and j indices (j < i for L (i,j) ( i,j ) and i < j for U (i,j) (u i,j )). Hence, we adopt the notation T (i,j) (t i,j ) to denote both L (i,j) ( i,j ) and U (i,j) (u i,j ), as we defined it earlier in equation (2.26).

Both iJDC +

LU-1 and iJDC + QR-1 algorithms consist of estimating the elementary triangular matrices U (i,j) (u i,j ). But the resulting parameters could be different. Additionally the procedures of estimating these matrices for both algorithms are identical.

Consequently, the principal problem is reduced to the estimation of two kinds of parameters, namely the elementary triangular matrices T (i,j) (t i,j ) and the Givens rotation matrices Q (i,j) (θ i,j ), by means of Jacobi-like optimization procedures.

Minimization with respect to

T (i,j) (t i,j )
In this section, we estimate T (i,j) (t i,j ) for i = j by solving the following optimization problem:

t i,j = argmin

t i,j Ψ + 1 ( BT (i,j) (t i,j )) (2.48)
The update of B is defined in equation (2.27) by B (new) = BT (i,j) (t i,j ). Inserting B (new) into the optimization (2.48), the optimal t i,j is the global minimum of the following local cost function:

Ψ + 1 (t i,j ) = K k=1 off [( B (new) ) 2 ] T C (k,-1) [( B (new) ) 2 ] 2 F (2.49)
2.3. The iJDC + LU-1 and iJDC + QR-1 algorithms Chapter 2

We define the K matrices C (k,new) ∈ R N ×N as follows:

C (k,new) = ( B (new) ) 2 T C (k,-1) ( B (new) ) 2 k ∈ {1, 2, . . . , K} (2.50)
The cost function (2.49) actually minimizes the total sum of squares of the off-diagonal components of C (k,new) . The explicit expression of ( B (new) ) 2 as a function of t i,j is provided in (2.29) in proposition 2.1. Inserting equation (2.29) into (2.50), we have:

C (k,new) = T (i,j) (t 2 i,j ) T C (k) T (i,j) (t 2 i,j ) 1 ○ + t i,j T (i,j) (t 2 i,j ) T c (k,1) e T j 2 ○ + t i,j e j c (k,2) T (i,j) (t 2 i,j ) 3 ○ + t 2 i,j c (k,3) e j e T j 4 ○ (2.51)
where

C (k) = A T C (k,-1) A is a (N ×N ) constant matrix, where c (k,1) = 2 A T C (k,-1) ( b i b j ) is a (N × 1) constant column vector, where c (k,2) = 2 ( b i b j ) T C (k,-1) A is a (1 × N )
constant row vector, and where

c (k,3) = 4 ( b i b j ) T C (k,-1) ( b i b j ) is a constant scalar.
The term 1 ○ in (2.51) transforms the j-th column vector and the j-th row vector of C (k) . The term 2 ○ in (2.51) is a zero matrix except its j-th column vector containing non-zero elements, while the term 3 ○ in (2.51) contains non-zero entries only on its j-th row vector. Eventually the term 4 ○ in (2.51) is a zero matrix except its (j, j)-th component being non-zero. In addition, C (k,new) = 1 ○+ 2 ○+ 3 ○+ 4 ○ is a (N ×N ) symmetric matrix. Hence (2.51) shows that only the j-th column vector and the j-th row vector of C (k,new) involve the parameter t i,j , while the other components remain constant. Therefore, considering the symmetry of C (k,new) , the minimization of the cost function (2.49) is equivalent to minimizing the total sum of the squares of the elements of the j-th column vectors of C (k,new) except their (j, j)-th elements with k ∈ {1, • • • , K}. The required components of C (k,new) can be expressed by the following proposition.

Proposition 2.2. The elements of the j-th column vector except the (j, j)-th entry of C (k,new) is a second degree polynomial function in t i,j as follows, for every n value different of j:

C (k,new) n,j = C (k) n,i t 2 i,j + c (k,1) n t i,j + C (k) n,j
(2.52)

where

C (k) n,i and C (k)
n,j are the (n, i)-th and (n, j)-th components of matrix C (k) , respectively, and c (k,1) n is the n-th element of vector c (k,1) .

The proof of this proposition is addressed in appendix A.2. Proposition 2.2 shows that the minimization of the cost function (2.49) can be expressed in the following compact matrix form:

Ψ + 1 (t i,j ) = K k=1 E (k) τ i,j 2 F = τ T i,j Q E τ i,j (2.53)
where Q E is defined as a (3 × 3) symmetric coefficient matrix:

Q E = K k=1 (E (k) ) T E (k) (2.54)
Chapter 2

2.3. The iJDC + LU-1 and iJDC + QR-1 algorithms with E (k) a ((N -1) × 3) matrix defined as follows: its first column contains the ith column of C (k) without the j-th element, the second column contains vector c (k,1) without the j-th entry and the third column contains the j-th column of C (k) without the j-th component. τ i,j is a 3-dimensional parameter vector defined as follows:

τ i,j = t 2 i,j , t i,j , 1 T (2.55)
Equation (2.53) shows that Ψ + 1 (t i,j ) is a fourth degree polynomial function. The global minimum t i,j can be obtained by computing the roots of its derivative and selecting the one yielding the smallest value of (2.53). Once the optimal t i,j is computed, B (new) is updated by means of equation (2.27) and the joint diagonalizer A is updated by computing ( B (new) ) 2 . The K matrices C (k,new) are computed by means of equation (2.51) and are used as the input matrices for the next Jacobi-like iteration. Then the same procedure is repeated to compute the next t i,j with another index couple (i, j). The processing of all the N (N -1) parameters t i,j with different (i, j) couples, more precisely i,j with i > j and u i,j with i < j, is called a LU sweep. The proposed iJDC + LU-1 algorithm requires several LU sweeps in practice.

Minimization with respect to

Q (i,j) (θ i,j )
Now we identify Q (i,j) (θ i,j ) with i < j and i, j ∈ {1, . . . , N } by solving the following optimization problem:

θ i,j = argmin θ i,j Ψ + 1 ( BQ (i,j) (θ i,j )) (2.56)
The update of B is now defined as follows:

B (new) = BQ (i,j) (θ i,j ) (2.57)
Similarly, putting (2.57) into the optimization (2.56), the optimal θ i,j is the global minimum of the following local cost function:

Ψ + 1 (θ i,j ) = K k=1 off [( B (new) ) 2 ] T C (k,-1) [( B (new) ) 2 ] 2 F (2.58)
We define the K matrices C (k,new) ∈ R N ×N as follows:

C (k,new) = ( B (new) ) 2 T C (k,-1) ( B (new) ) 2 , k ∈ {1, 2, . . . , K} (2.59)
The cost function (2.58) minimizes the total sum of squares of the off-diagonal components of C (k,new) . The explicit expression of ( B (new) ) 2 as a function of θ i,j now can be rewritten as shown in the following proposition.

Proposition 2.3. A (new) = ( B (new) ) 2 = ( BQ (i,j) (θ i,j )) 2 can be written as a function of θ i,j as follows:

A (new) = ( B (new) ) 2 = B 2 (Q (i,j) (θ i,j )) 2 + sin(2θ i,j )( b i b j )(e T i -e T j ) (2.60)
where b i and b j denote the i-th and j-th column vectors of B, respectively, and where e i and e j are the i-th and j-th column vectors of the identity matrix I N , respectively.

The iJDC +

LU-1 and iJDC + QR-1 algorithms Chapter 2

The proof of proposition 2.3 is provided in appendix A.3. Inserting (2.60) into (2.59), we obtain:

C (k,new) =[(Q (i,j) (θ i,j )) 2 ] T C (k) (Q (i,j) (θ i,j )) 2 1 ○ + sin(2θ i,j )[(Q (i,j) (θ i,j )) 2 ] T c (k,1) (e T i -e T j ) 2 ○ + sin(2θ i,j )(e i -e j ) c (k,2) (Q (i,j) (θ i,j )) 2 3 ○ + sin 2 (2θ i,j ) c (k,3) (e i -e j )(e T i -e T j ) 4 ○ (2.61) where C (k) = A T C (k,-1) A is a (N ×N ) constant matrix, where c (k,1) = A T C (k,-1) ( b i b j ) is a (N × 1) constant column vector, where c (k,2) = ( b i b j ) T C (k,-1) A is a (1 × N ) constant row vector, and where c (k,3) = ( b i b j ) T C (k,-1) ( b i b j ) is a constant scalar. The term 1
○ in (2.61) transforms the i-th and j-th column vectors and the i-th and j-th row vectors of C (k) . The term 2 ○ in (2.61) is a zero matrix except its i-th and j-th column vectors containing non-zero elements, while the term 3 ○ in (2.61) contains non-zero entries only on its i-th and j-th row vectors. And the term 4 ○ in (2.61) is a zero matrix except its (i, i)-th, (j, j)-th, (i, j)-th and (j, i)-th components being non-zero.

C (k,new) = 1 ○+ 2 ○+ 3 ○+ 4
○ is a (N ×N ) symmetric matrix. Hence (2.61) shows that only the i-th and j-th column vectors and the i-th and j-th row vectors of C (k,new) involve the parameter θ i,j , while the other components remain constant. It is noteworthy that the (i, j)-th and (j, i)-th components are twice affected by the transformation. Considering the symmetry of C (k,new) , we propose to minimize the sum of the squares of the (i, j)-th entries of the K matrices C (k,new) , instead of minimizing all the off-diagonal entries. Although minimizing this quantity is not equivalent to minimizing the original cost function (2.58), such a simplified minimization scheme is commonly adopted by many algorithms, such as in [START_REF] Wang | Nonorthogonal joint diagonalization algorithm based on trigonometric parameterization[END_REF], Souloumiac, 2009, Guo et al., 2010, Luciani and Albera, 2014]. Therefore, u i,j is estimated by minimizing an approximate of cost function (2.58), denoted by Ψ + 1 (θ i,j ). The (i, j)-th component of C (k,new) is expressed in the following proposition. Proposition 2.4. The (i, j)-th entry of C (k,new) can be expressed as a function of θ i,j as follows:

C (k,new) i,j = -sin 2 (2θ i,j ) c (k,3) + sin 2 (θ i,j )( C (k) i,i cos 2 (θ i,j ) + C (k) j,i sin 2 (θ i,j )) + cos 2 (θ i,j )( C (k) i,j cos 2 (θ i,j ) + C (k) j,j sin 2 (θ i,j )) + sin(2θ i,j )( c (k,1) i cos 2 (θ i,j ) + c (k,1) j sin 2 (θ i,j )) -sin(2θ i,j )( c (k,2) j cos 2 (θ i,j ) + c (k,2) i sin 2 (θ i,j )) (2.62) where C (k) i,i , C (k) j,j , C (k)
i,j and C (k) j,i are the (i, i)-th, (j, j)-th, (i, j)-th and (j, i)-th components of matrix C (k) , respectively, and where c (k,q) i and c

(k,q) j are the i-th and j-th elements of vector c (k,q) with q ∈ {1, 2}, respectively.

The proof of proposition 2.4 is given in appendix A.4. In order to simplify the notation of (2.62), we resort to the Weierstrass change of variable: ϑ i,j = tan(θ i,j ). Then we obtain:

sin(2θ i,j ) = 2ϑ i,j 1 + ϑ 2 i,j , cos(2θ i,j ) = 1 -ϑ 2 i,j 1 + ϑ 2 i,j , sin 2 (θ i,j ) = ϑ 2 i,j 1 + ϑ 2 i,j , cos 2 (θ i,j ) = 1 1 + ϑ 2 i,j
(2.63) Chapter 2

The iJDC +

LU-1 and iJDC + QR-1 algorithms By substituting (2.63) into (2.62), we obtain an alternative expression of the (i, j)-th entry of C (k,new) which is described in the following proposition.

Proposition 2.5. The (i, j)-th entry of C (k,new) can be expressed by a rational function of ϑ i,j as follows:

C (k,new) i,j = f (k) 4 ϑ 4 i,j + f (k) 3 ϑ 3 i,j + f (k) 2 ϑ 2 i,j + f (k) 1 ϑ i,j + f (k) 0 (1 + ϑ 2 i,j ) 2
(2.64)

where f

(k) 4 = C (k) j,i , f (k) 3 = -2 c (k,1) i , f (k) 2 = C (k) i,i + C (k) j,j + 2 c (k,2) j -4 c (k,3) , f (k) 1 = 2 c (k,2) i - c (k,1) j and f (k) 0 = C (k)
j,j . The proof of proposition 2.5 is quite easy and therefore omitted. By means of proposition 2.5, the approximate cost function

Ψ + 1 (θ i,j ) is then transformed into Ψ + 1 (ϑ i,j ). Recall that Ψ + 1 (ϑ i,j
) is the total sum of the squares of the (i, j)-th entries of the K matrices C (k,new) . Moreover equation (2.64) easily shows that Ψ + 1 (ϑ i,j ) is a rational function in ϑ i,j , where the degrees of the numerator and the denominator are 8 and 8, respectively. Ψ + 1 (ϑ i,j ) can be expressed in the following compact matrix form:

Ψ + 1 (ϑ i,j ) = K k=1 (f (k) ) T ϑ i,j 2 F = ϑ T i,j Q F ϑ i,j (2.65)
where Q F is a (5 × 5) symmetric coefficient matrix defined as follows:

Q F = K k=1 f (k) (f (k) ) T (2.66)
where f (k) is a 5-dimensional vector:

f (k) = f (k) 4 , f (k) 3 , f (k) 2 , f (k) 1 , f (k) 0 T (2.67)
and where ϑ i,j is a 5-dimensional parameter vector defined as follows:

ϑ i,j = 1 (1 + ϑ 2 i,j ) 2 ϑ 4 i,j , ϑ 3 i,j , ϑ 2 i,j , ϑ i,j , 1 T (2.68)
The global minimum ϑ i,j can be obtained by computing the roots of the derivative of the cost function(2.65) and selecting the one yielding the smallest value of (2.65). Once ϑ i,j is obtained, θ i,j can be computed from the inverse tangent function θ i,j = arctan(ϑ i,j ). It is noteworthy that the found θ i,j cannot guarantee to decrease the original cost function (2.58). If θ i,j leads to an increase of (2.58), we reset θ i,j = 0, thus the value of the original cost function remains unchanged. Otherwise, B is updated as described in (2.57) and the joint diagonalizer A is updated by computing ( B (new) ) 2 . The same procedure will be repeated to compute θ i,j with the next (i, j) couple. The processing of all the N (N-1)/2 parameters θ i,j and also the other N (N -1)/2 parameters u i,j is called a QR sweep. Several QR sweeps yield the proposed iJDC + QR-1 algorithm.

Both the iJDC + LU-1 and iJDC + QR-1 algorithms can be stopped when the value of the global cost function (2.47) or its relative change between two successive sweeps fall below a fixed small positive threshold. Such a stopping criterion is guaranteed to be met since the cost function is non-increasing in each Jacobi-like sweep.

The iJDC +

LU-1 and iJDC + QR-1 algorithms Chapter 2

Practical issues

Nonnegativity relaxation

As we discussed in section 2.2.3, if each input matrix C (k) is almost exactly jointly diagonalizable due to a high SNR value, the explicit nonnegativity constraint by means of a square change of variable could be unnecessary. Consequently, the nonnegativity constraint can be relaxed by directly decomposing A using the LU factorization in equation (2.43) or the QR factorization as follows:

A = N i=1 N j=i+1 Q (i,j) (θ i,j ) N i=1 N j=i+1 U (i,j) (u i,j ) def = i<j Q (i,j) (θ i,j ) i<j T (i,j) (t i,j ) (2.69)
respectively. By means of the JDC criterion Ψ 1 (A) in (2.46), each T (i,j) (t i,j ) and Q (i,j) (θ i,j ) can be estimated by solving the two following Jacobi-like optimization problems, respectively:

t i,j = argmin t i,j Ψ 1 ( AT (i,j) (t i,j )) (2.70) θ i,j = argmin θ i,j Ψ 1 ( AQ (i,j) (θ i,j )) (2.71)
where A is the estimate of A at the (it -1)-th iteration. It is easy to show that the right hand side of (2.70) is a second degree polynomial function in t i,j . Then the global optimal t i,j can be directly expressed as follows:

t i,j = - k n =j C (k) i,n C (k) j,n k n =j C (k) i,n C (k) i,n
(2.72)

where C (k) = A T C (k,-1) A is a (N × N ) constant matrix, and where C (k)

i,n and C

(k) j,n are the (i, n)-th and (j, n)-th components of C (k) , respectively. In addition, the right hand side of (2.71) can be written as the following quadratic form:

Ψ 1 (θ i,j ) = θ T i,j (G T G)θ i,j (2.73)
where θ i,j is a 2-dimensional parameter vector defined by θ i,j = [cos (2θ i,j ), sin (2θ i,j )] T , and where G ∈ R K×2 is a coefficient matrix, whose (k, 1)-th and (k, 2)-th components are equal to

C (k) j,j -C (k) i,i and C (k) i,j + C (k)
j,i , respectively. The vector θ i,j can be chosen as the unit-norm eigenvector corresponding to the smaller eigenvalue of the matrix G T G. Then θ i,j can be computed from θ i,j accordingly. The derivations of equations (2.72) and (2.73) are straightforward and therefore omitted.

In fact, the ways of estimating the two sets of parameters { i,j , u i,j } in equation (2.43), and {θ i,j , u i,j } in equation (2.69) are identical to those of Afsari's LUJ1D and QRJ1D methods [Afsari, 2006] performing on the inverted input matrices C (k,-1) , respectively. In practice, for iJDC + LU-1 and iJDC + QR-1 , in each Jacobi-like iteration for computing an elementary triangular matrix, we suggest to compute t i,j by equation (2.72) first. If all the elements in the j-th column vector of AT (i,j) (t i,j ) have the same sign ς, the j-th column vector of A (new) is replaced by ς( a j + t i,j a i ), while the other column vectors remain Chapter 2

The iJDC +

LU-1 and iJDC + QR-1 algorithms unchanged, where a i and a j are the i-th and j-th column vectors of A, respectively. Otherwise, t i,j is computed by minimizing (2.53). Then A is updated by means of equation (2.29). Regarding the Givens rotation matrices Q (i,j) (θ i,j ) in the iJDC + QR-1 algorithm, θ i,j is computed in the same manner by minimizing either (2.73) or (2.65).

Then A (new) is computed accordingly.

Row balancing

Afsari reported in [Afsari, 2006] that if the row vectors of matrices

C (k) (k ∈ {1, • • • , K})
are not balanced in their norms, the computation of the parameters could be inaccurate.

In order to cope with this effect, we apply Afsari's row balancing scheme every few sweeps. Such a scheme updates each

C (k) by C (k,new) = Λ C (k) Λ and A by A (new) = AΛ using a diagonal matrix Λ ∈ R N ×N +
, whose diagonal elements are defined as follows:

Λ n,n = 1 K k=1 C (k) n,: 2 , n ∈ {1, 2, • • • , N } (2.74)
where

C (k)
n,: denotes the n-th row vector of matrix C (k) .

Discussions

The effectiveness of the iJDC + LU-1 and iJDC + QR-1 algorithms relies on the two assumptions that A is a nonsingular square matrix and that the main diagonals of {D (k) } do not contain any zero entry. In ICA, the first nonsingularity assumption of the mixing matrix A is mostly satisfied. In the following we verify the validity of the second assumption.

Provided that all the sources are non-Gaussian, which is often the case in practice [Comon and Jutten, 2010], we can resort to use the FO cumulant array of the observation data. As we mentioned in section 1.3.1, the jointly diagonalizable matrices {C (k) } can be built by stacking the matrix slices of the cumulant array. Then the (N 2 × N ) matrix D whose k-th row vector contains the diagonal components of D (k) , can be expressed by D

= (A A)C (4) s , where C (4) s = diag [C 1,1,1,1,s , . . . , C N,N,N,N,s ] is a (N × N ) diagonal matrix with C n,
n,n,n,s being the FO cumulant of the n-th source, n ∈ {1, . . . , N }, and where denotes the Khatri-Rao product. Suppose that A is a nonsingular dense mixing matrix, then D is also nonsingular and contains non-zeros components. A dense mixing matrix means that the sources are highly mixed, therefore the components of A have relative large values. A dense mixing matrix is often encountered in BSS problems, such as in hyperspectral unmixing [Miao andQi, 2007, Chan et al., 2009] and in separating MRS signals [Moussaoui, 2005]. Consequently, the nonzero assumption is reasonable for some BSS applications.

One drawback brought by these assumptions is that the mixing matrix (and the joint transformation matrix) A is constrained to be square, since a non-square matrix A ∈ R N ×P + with N > P will lead to the jointly diagonalizable matrices C (k) = AD (k) A T being rank deficient and non-invertible. In ICA, when such a non-square matrix A is encountered, we can compress A by means of a nonnegative matrix W ∈ R P ×N + such that the resulting matrix A + = WA is a nonnegative square matrix. W can be computed by means of the nonnegative compression algorithm that we proposed in [START_REF] Wang | Nonnegative compression for semi-nonnegative independent component analysis[END_REF], which transforms the classical prewhitening matrix [Belouchrani et al., 1997] 2.4. The JDC + LU-1 and JDC + LU-2 algorithms Chapter 2 into a nonnegative one by linear mappings. Such a method is described in appendix B.

Then the iJDC + LU-1 and iJDC + QR-1 algorithms are used to compute the compressed mixing matrix A + .

Finally, the pseudo-codes for the iJDC + LU-1 and iJDC + QR-1 algorithms are given in the appendix C.3 and C.4, respectively.

The JDC +

LU-1 and JDC + LU-2 algorithms

Algorithm derivation

In some practical JDC problems, the main diagonals of matrices D (k) could contain zero entries, which breaks the second assumption of the iJDC + LU-1 and iJDC + QR-1 algorithms. Therefore such an assumption could limit the applicability of iJDC + LU-1 and iJDC + QR-1 . Moreover, even though this assumption is satisfied, in some cases, some jointly diagonalizable matrices C (k) may be ill-conditioned. The inverse preprocessing of these matrices could generate bias in the following JDC procedure. Consequently, the purpose of this section is to develop Jacobi-like algorithms that require less assumptions. More precisely, we allow the diagonal matrices D (k) contain zero components on their main diagonals and discard the preprocessing step of inverting each C (k) . We still assume that the nonnegative joint transformation matrix A is nonsingular. In addition, we would like to study the performance of the Jacobi-like procedure on difference cost functions. Hence, by considering both the indirect-fit JDC criteria (1.27) and (1.28) and the LU matrix factorization, we derive two algorithms, namely JDC + LU-1 and JDC + LU-2 . By inserting the square change of variable A = B 2 into the indirect-fit JDC criteria (1.27) and (1.28), we obtain the following two nonnegativity constrained JDC criteria, respectively:

Ψ + 1 (B) = K k=1 off{(B 2 ) -1 C (k) (B 2 ) -T } 2 F
(2.75)

Ψ + 2 (B) = K k=1 C (k) -B 2 diag{(B 2 ) -1 C (k) (B 2 ) -T }(B 2 ) T 2 F (2.76)
By means of the LU matrix factorization of B and the notation B = i =j T (i,j) (t i,j ) defined in equation (2.26), the minimization of (2.75) and (2.76) with respect to B is converted to the estimation of the set of N (N -1) parameters t i,j . Therefore, each parameter t i,j can be estimated sequentially using one of the two following Jacobi-like optimization procedures:

t i,j = argmin t i,j Ψ + 1 (B (it-1) T (i,j) (t i,j )) (2.77) t i,j = argmin t i,j Ψ + 2 (B (it-1) T (i,j) (t i,j )) (2.78)
As in the previous sections, let A and B denote A (it-1) and B (it-1) , respectively, before estimating the parameter t i,j . Let A (new) and B (new) stand for A and B updated by T (i,j) (t i,j ), respectively. The update of B is defined by B (new) = BT (i,j) (t i,j ). Inserting B (new) into the Jacobi-like optimizations (2.77) and (2.78), the optimal t i,j can be chosen Chapter 2

2.4. The JDC + LU-1 and JDC + LU-2 algorithms as the global minimum of the following two local cost functions, respectively:

Ψ + 1 (t i,j ) = K k=1 off ( B (new) ) 2 -1 C (k) ( B (new) ) 2 -T 2 F (2.79) Ψ + 2 (t i,j ) = K k=1 C (k) -( B (new) ) 2 × diag ( B (new) ) 2 -1 C (k) ( B (new) ) 2 -T ( B (new) ) 2 T 2 F (2.80)
We define the K matrices C (k,new) ∈ R N ×N as follows:

C (k,new) = ( B (new) ) 2 -1 C (k) ( B (new) ) 2 -T , k ∈ {1, 2, . . . , K} (2.81)
According to equations (2.79), (2.80) and (2.81), the computation of the cost functions Ψ + 1 (t i,j ) and Ψ + 2 (t i,j ) requires to express the components of the K matrices C (k,new) as an explicit function of t i,j . Therefore, it is necessary to write [( B (new) ) 2 ] -1 as a function of t i,j . The explicit expression of ( B (new) ) 2 as a function of t i,j is provided in equation (2.29) in proposition 2.1. In the right hand side of (2.29), the first term of the sum is a nonsingular matrix and the second term is a rank-1 matrix. The inverse of the sum of such two matrices can be computed by the Sherman-Morrison formula [Sherman andMorrison, 1950, Bartlett, 1951]:

Theorem 2.1. Suppose that R ∈ R N ×N is a nonsingular square matrix and u, v ∈ R N are two column vectors satisfying 1 + v T R -1 u = 0, then: (R + uv T ) -1 = R -1 - R -1 uv T R -1 1 + v T R -1 u (2.82)
Suppose that matrix B, the two vectors 2 t i,j ( b i b j ) and e j satisfy the conditions of the Sherman-Morrison formula, the expression of [( B (new) ) 2 ] -1 as a function of the parameter t i,j has the following form:

( B (new) ) 2 -1 = T (i,j) (-t 2 i,j ) Q( B 2 ) -1 (2.83)
where Q ∈ R N ×N is defined as follows:

Q = I N - 2 t i,j 1 + 2β j t i,j βe T j (2.84)
where I N ∈ R N ×N is the identity matrix, where β ∈ R N is a column vector defined as follows:

β = ( B 2 ) -1 ( b i b j ) (2.85)
and where β j is the j-th element of vector β. Inserting equation (2.83) into (2.81), C (k,new) can be rewritten as follows:

C (k,new) = T (i,j) (-t 2 i,j ) Q C (k) Q T T (i,j) (-t 2 i,j ) T (2.86)
where

C (k) = A -1 C (k) A -T is a (N × N ) constant matrix.
Then through a straightforward computation of (2.86), each component of C (k,new) can be expressed as a function of the parameter t i,j as described in the following two propositions.

2.4. The JDC + LU-1 and JDC + LU-2 algorithms Chapter 2

Proposition 2.6. Each non-diagonal component of C (k,new) is a rational function in t i,j :

C (k,new) m,n m =n = E (k,3) m,n t 3 i,j + E (k,2) m,n t 2 i,j + E (k,1) m,n t i,j + E (k,0) m,n (1 + 2β j t i,j ) 2 (2.87)
where

E (k,3) m,n , E (k,2) m,n , E (k,1) m,n and E (k,0)
m,n are the (m, n)-th components of the (N × N ) symmetric coefficient matrices E (k,3) , E (k,2) , E (k,1) and E (k,0) , respectively. These coefficients are defined as follows:

E (k,3) m,n =        2( C (k) j,j β m -C (k) m,j β j ) if n = i, 1 ≤ m = i ≤ N E (k,3) n,m if m = i, 1 ≤ n = i ≤ N 0 otherwise
(2.88)

E (k,2) m,n =                        4( C (k) j,j β m β n + C (k) m,n β 2 j -( C (k) m,j β n + C (k) j,n β m )β j ) if 1 ≤ m < n ≤ N, m = i, n = i 4( C (k) j,j β m β n + C (k) m,n β 2 j -( C (k) m,j β n + C (k) j,n β m )β j ) -C (k) m,j if n = i, 1 ≤ m < i E (k,2) n,m if 1 ≤ n < m ≤ N 0 otherwise
(2.89)

E (k,1) m,n =        4 C (k) m,n β j -2( C (k) m,j β n + C (k) j,n β m ) if 1 ≤ m < n ≤ N E (k,1) n,m if 1 ≤ n < m ≤ N 0 otherwise
(2.90)

E (k,0) m,n = C (k) m,n if 1 ≤ m = n ≤ N 0 otherwise (2.91)
where C

(k) m,n is the (m, n)-th component of the matrix C (k) , and where β m is the m-th element of the vector β which is defined in equation (2.85).

Proposition 2.7. Let us define K column vectors d (k) ∈ R N with k ∈ {1, . . . , K}, each containing the diagonal elements of C (k,new) , then the n-th element of d (k) , denoted by d

(k)
n , can be expressed as a rational function in t i,j as follows:

d (k) n = f (k,4) n t 4 i,j + f (k,3) n t 3 i,j + f (k,2) n t 2 i,j + f (k,1) n t i,j + f (k,0) n (1 + 2β j t i,j ) 2 (2.92)
where f

(k,4) n , f (k,3) n , f (k,2) n , f (k,1) n and f (k,0) n
are the n-th elements of the N -dimensional coefficient vectors f (k,4) , f (k,3) , f (k,2) , f (k,1) and f (k,0) , respectively, with n ∈ {1, . . . , N }. These coefficients are defined as follows:

f (k,4) n = C (k) j,j if n = i 0 otherwise (2.93) f (k,3) n = 4 C (k) j,j β i -2( C (k) i,j + C (k) j,i )β j if n = i 0 otherwise (2.94)
Chapter 2

2.4. The JDC + LU-1 and JDC + LU-2 algorithms

f (k,2) n =                    4( C (k) n,n β 2 j -( C (k) n,j + C (k) j,n )β n β j + C (k) j,j β 2 n ) if n = i, n = j 4 C (k) i,i β 2 j -( C (k) i,j + C (k) j,i )(1 + 4β i β j ) + 4 C (k) j,j β 2 i if n = i 0 otherwise
(2.95)

f (k,1) n = 4 C (k) n,n β j -2( C (k) n,j + C (k) j,n )β n if n = j 0 otherwise (2.96) f (k,0) n = C (k) n,n
(2.97)

where C

(k) i,j is the (i, j)-th component of the matrix C (k) , and where β i is the i-th element of the vector β which is defined in equation (2.85).

The proofs of propositions 2.6 and 2.7 are provided in the appendix A.5. Then we can derive the proposed JDC + LU-1 and JDC + LU-2 algorithms.

The JDC +

LU-1 algorithm The cost function Ψ + 1 (t i,j ) defined in (2.79) is equal to the total sum of the squares of the non-diagonal components of C (k,new) defined in proposition 2.6, therefore it can be expressed in the following compact matrix form:

Ψ + 1 (t i,j ) = K k=1 E (k) τ i,j 2 F = τ T i,j Q E τ i,j (2.98)
where Q E ∈ R 4×4 is a symmetric coefficient matrix defined by:

Q E = K k=1 (E (k) ) T E (k) (2.99)
where E (k) is a (N 2 × 4) matrix defined as follows:

E (k) = vec (E (k,3) ), vec (E (k,2) ), vec (E (k,1) ), vec (E (k,0) ) (2.100)
and where τ i,j is defined as the following 4-dimensional parameter vector :

τ i,j = 1 (1 + 2β j t i,j ) 2 [t 3 i,j , t 2 i,j , t i,j , 1] T (2.101) Equation (2.98) shows that Ψ + 1 (t i,j
) is a rational function, where the degrees of the numerator and the denominator are 6 and 4, respectively. The global minimum t i,j can be obtained by computing the roots of the numerator of its derivative,which is a polynomial of degree 9, and selecting the one yielding the smallest value of (2.98). In practice, the roots of a polynomial of degree 9 can be obtained by computing the eigenvalues of the corresponding companion matrix of dimension (9 × 9) [Golub and Van Loan, 1996].

The JDC +

LU-1 and JDC + LU-2 algorithms Chapter 2 2.4.1.2 The JDC + LU-2 algorithm According to proposition 2.7, d (k) can be expressed in the following compact matrix form:

d (k) = F (k) ζ i,j (2.102)
where F (k) ∈ R N ×5 is a coefficient matrix:

F (k) = f (k,4) , f (k,3) , f (k,2) , f (k,1) , f (k,0) (2.103)
and where ζ i,j ∈ R 5 is a parameter vector defined as follows:

ζ i,j = 1 (1 + 2β j t i,j ) 2 [t 4 i,j , t 3 i,j , t 2 i,j , t i,j , 1] T (2.104)
The use of vec (.) operator in criterion (2.80) leads to the following equivalent form of Ψ + 2 (t i,j ):

Ψ + 2 (t i,j ) = K k=1 c (k) -( B (new) ) 2 ( B (new) ) 2 d (k) 2 F (2.105)
where c (k) = vec (C (k) ), and where denotes the Khatri-Rao product. Let us consider the following proposition:

Proposition 2.8. According to equation (2.29) and the definition of Khatri-Rao product, the following equality holds:

( B (new) ) 2 ( B (new) ) 2 = ( A A) + G ζ i,j e T j (2.106)
where G is a (N 2 ×5) matrix, whose column vectors are defined as follows:

g 1 = vec ( a i a T i ) (2.107) g 2 = vec (2 a i ( b i b j ) T + 2 ( b i b j ) a T i )
(2.108)

g 3 = vec ( a i a T j + a j a T i + 4 ( b i b j )( b i b j ) T ) (2.109) g 4 = vec (2 a j ( b i b j ) T + 2 ( b i b j ) a T j )
(2.110)

g 5 = 0 (2.111)
and where a i and a j are the i-th and j-th column vectors of A, respectively.

The proof of proposition 2.8 is given in appendix A.6. By inserting equations (2.102) and (2.106) into (2.105), the cost function Ψ + 2 (t i,j ) can be rewritten as follows:

Ψ + 2 (t i,j ) = K k=1 c (k) -G (k) A ζ i,j 2
F (2.112)
where

G (k)
A is a (N 2 × 5) matrix defined as follows:

G (k) A = ( A A)F (k) + f (k,0) j G (2.113)
Now, let us define a (N 2 × 5) matrix H (k) as follows:

H (k) = c (k) (4 β 2 j e T 3 + 4 β j e T 4 + e T 5 ) -G (k) A (2.114)
Chapter 2

2.4. The JDC + LU-1 and JDC + LU-2 algorithms where e n , n ∈ {3, 4, 5} is the n-th column vector of the (5 × 5) identity matrix. Then the criterion Ψ + 2 (t i,j ) can be expressed in a compact matrix form as follows: 2.115) where:

Ψ + 2 (t i,j ) = K k=1 H (k) ζ i,j 2 F = ζ T i,j Q H ζ i,j (
Q H = K k=1 (H (k) ) T H (k) (2.116)
is a (5 × 5) symmetric coefficient matrix. Equation (2.115) shows that the cost function Ψ + 2 (t i,j ) is a rational function, where the degrees of the numerator and the denominator are 8 and 4, respectively. The global minimum t i,j can be obtained by computing the roots of its derivative and selecting the one yielding the smallest value of (2.115).

Once the optimal t i,j is computed by either JDC + LU-1 or JDC + LU-2 , B (new) and A (new) are computed as described in proposition 2.1. Then the Jacobi-like procedure continues to compute t i,j with the next (i, j) index couple. The processing of all the N (N -1) parameters is called a sweep. In practice, it is necessary to perform several sweeps in order to ensure the convergence of this sequential optimization procedure. Such procedures can be terminated when the values of (2.98) and (2.115), or their relative decreases between two successive sweeps fall below a specified small positive threshold. Such a stopping criterion is guaranteed to be met since the criterions (2.98) and (2.115) are guaranteed to decrease (or at least to not increase) after each Jacobi-like sweep.

Generalization to the non-square case

In practice a JDC problem often involves a non-square joint transformation matrix A ∈ R N ×P + and K diagonal matrices D (k) ∈ R P ×P with N > P . In order to compute a nonsquare full column rank matrix A by minimizing the cost functions (1.27) and (1.28), the inverse operator in (1.27) and (1.28) should be replaced by the pseudo-inverse operator. However, the pseudo-inverse is more difficult to handle than the inverse operation. In order to avoid computing the pseudo-inverse, let us firstly define the following two (N ×N ) matrices:

Å = [A, e P +1 , . . . , e N ] = [a 1 , a 2 , . . . , a P , e P +1 , . . . , e N ]
(2.117)

D(k) = D (k) 0 P,(N -P ) 0 (N -P ),P 0 (N -P ),(N -P ) (2.118)
where a p is the p-th column vector of A, where e p is the p-th column vector of the identity matrix I N , and where 0 Y,Z denotes a (Y × Z) zero matrix. Then we obtain the following equality: (2.121)

C (k) = A D (k) A T = Å D(k) ÅT (2.
Ů = U 0 P ×(N -P ) 0 (N -P )×P I N -P (2.122)
where B ∈ R N ×P , where L ∈ R N ×P is a unit lower triangular matrix, where U ∈ R P ×P is a unit upper triangular matrix, and where I N -P ∈ R (N -P )×(N -P ) is an identity matrix. In addition B, L and U satisfy B = LU . If A = B 2 , it is easy to prove that Å = B 2 = ( L Ů ) 2 . Furthermore, L and Ů admit the following elementary triangular factorization, respectively:

L = P j=1 N i=j+1 L(i,j) ( i,j ), Ů = P i=1 P j=i+1 Ů (i,j) (u i,j ) (2.123)
where L(i,j) ( i,j ) ∈ R N ×N and Ů (i,j) (u i,j ) ∈ R N ×N are elementary lower and upper triangular matrices, respectively. Therefore B can be directly estimated by the proposed JDC + LU-1 and JDC + LU-2 methods. The total number of free parameters of B is equal to

P (N -1).
Then the final output matrix B ∈ R N ×P can be obtained be choosing the first P column vectors of B ∈ R N ×N .

Practical issues

As we mentioned in the previous sections, when the model error of each jointly diagonalizable matrix C (k) is small, the explicit nonnegativity constraint on the matrix A could be needless. We can directly decompose A using the LU factorization as presented in equation (2.43). By means of the JDC criteria Ψ 1 (A) (1.27) and Ψ 2 (A) (1.28), each T (i,j) (t i,j ) in decomposition (2.43) can now be estimated by solving one of the two following Jacobi-like optimization problems, respectively:

t i,j = argmin t i,j Ψ 1 ( AT (i,j) (t i,j )) (2.124) t i,j = argmin t i,j Ψ 2 ( AT (i,j) (t i,j )) (2.125)
where A is the estimate of A at the (it-1)-th iteration. By a straightforward computation of (2.124), its global optimal t i,j is the minimum of a second degree polynomial, and can be expressed as follows:

t i,j = k n =i C (k) i,n C (k) j,n k n =i C (k) j,n C (k) j,n
(2.126)

where are the (i, n)-th and (j, n)-th components of C (k) , respectively. Regarding (2.125), the function to be minimized can be approximated by a fourth degree polynomial function Chapter 2

C (k) = A -1 C (k) A -T is a (N × N ) constant matrix
2.5. The JDC + ADMM-3 algorithm [Afsari, 2006]:

Ψ 2 ( AT (i,j) (t i,j )) ≈ 4 K k=1 ( C (k) j,j ) 2 t 4 i,j -8 K k=1 C (k) j,j C (k) i,j t 3 i,j + K k=1 [4( C (k) i,j ) 2 + 2( C (k) j,j ) 2 ] t 2 i,j -4 K k=1 C (k) j,j C (k) i,j t i,j + K k=1 ( C (k) i,j ) 2
(2.127) where C

(k) i,j and C

(k) j,j are the (i, j)-th and (j, j)-th components of the constant matrix C (k) = A -1 C (k) A -T , respectively. The global minimum t i,j of equation ( 2.127) can be obtained by computing the roots of its derivative. The derivations of equations (2.126) and (2.127) are straightforward therefore omitted. As a matter of fact, the estimations of t i,j by using equations (2.126) and (2.127) are identical to the Afsari's LUJ1D and LUJ2D algorithms [Afsari, 2006], respectively, with a change of sign in t i,j . Practically, in each Jacobi-like iteration of JDC + LU-1 , we suggest to compute t i,j by equation (2.126) first. If all the elements in the j-th column vector of AT (i,j) (t i,j ) have the same sign ς, the j-th column vector of A (new) is replaced by ς( a j + t i,j a i ), while the other column vectors remain unchanged, where a i and a j are the i-th and j-th column vectors of A, respectively. Otherwise, t i,j is computed by minimizing (2.98). Then A is updated by means of equation (2.29). Regarding the JDC + LU-2 algorithm, t i,j is computed in the same manner by minimizing either (2.127) or (2.115). Then A is computed correspondingly.

We also utilize the row balancing scheme for the proposed JDC + LU-1 and JDC + LU-2 algorithms, as described in equation (2.74), where the matrix

C (k) is replaced by C (k) = A -1 C (k) A -T .
For conclusion, the pseudo-codes for the JDC + LU-1 and JDC + LU-2 algorithms are presented in the appendix C.5 and C.6, respectively.

The JDC + ADMM-3 algorithm

The previous three sections presented five Jacobi-like NJDC algorithms. The purpose of this section is to introduce another algorithm based on the Alternating Direction Method of Multipliers (ADMM) approach for solving the NJDC model.

Algorithm derivation

We propose to impose the nonnegativity constraint in the direct-fit JDC criterion (1.29) as follows:

Ψ + 3 (A, {D (k) }) = 1 2 K k=1 C (k) -AD (k) A T 2 F (2.128) subject to A ≥ 0
In order to facilitate an efficient use of alternating minimization, we introduce three

auxiliary variable A 1 ∈ R N ×P , A 2 ∈ R N ×P and U ∈ R N ×P .
Then the NJDC problem can be solved by minimizing the following cost function:

Ψ + 3 (A 1 , A 2 , {D (k) }) = 1 2 K k=1 C (k) -A 1 D (k) A T 2 2 F (2.129) subject to A 1 -U = 0, A 2 -U = 0, U ≥ 0 2.5. The JDC + ADMM-3 algorithm Chapter 2
where {D (k) } ⊂ R P ×P are a set of K diagonal matrices. The augmented Lagrangian function of (2.129) is expressed as follows:

L A (A 1 , A 2 , U , {D (k) }, Π 1 , Π 2 ) = Ψ + 3 (A 1 , A 2 , {D (k) }) + f (U ) + Π 1 • (A 1 -U ) + Π 2 • (A 2 -U ) + α 2 ||A 1 -U || 2 F + β 2 ||A 2 -U || 2 F (2.130)
where f (•) is the indicator function on the set of (N × P ) nonnegative matrices, where Π 1 ∈ R N ×P and Π 2 ∈ R N ×P are the Lagrangian multipliers, and where α > 0 and β > 0 are penalty parameters for the constraints A 1 -U = 0 and A 2 -U = 0, respectively. For the matrices Y and Z of the same dimension, the scalar product • is defined as the sum of all the element-wise products:

Y • Z = i j y i,j z i,j (2.131)
where y i,j and z i,j are the (i, j)-th components of Y and Z, respectively. The ADMM for (2.129) is derived by successively minimizing the augmented Lagrangian function (2.130) with respect to A 1 , A 2 , U and {D (k,it) }, one at a time while fixing others at their most recent values, for example:

A (it+1) 1 = argmin A 1 L A (A 1 , A (it) 2 , U (it) , {D (k,it) }, Π (it) 1 , Π (it) 2 )
(2.132)

A (it+1) 2 = argmin A 2 L A (A (it+1) 1 , A 2 , U (it) , {D (k,it) }, Π (it) 1 , Π (it)
2 ) (2.133)

U (it+1) = argmin U L A (A (it+1) 1 , A (it+1) 2 , U , {D (k,it) }, Π (it)
1 , Π

2 ) (2.134)

{D (k,it+1) } = argmin D (k) L A (A (it+1) 1 , A (it+1) 2 , U (it+1) , {D (k) }, Π (it) 1 , Π (it)
2 ) (2.135) and then updating the multipliers Π 1 and Π 2 , where it denotes the iteration number. Specifically, these steps can be written in closed form as follows:

A (it+1) 1 = K k=1 C (k) A (it) 2 D (k,it) + αU (it) -Π (it) 1 × K k=1 D (k,it) (A (it) 2 ) T A (it) 2 D (k,it) + αI P -1
(2.136)

A (it+1) 2 = K k=1 C (k) A (it+1) 1 D (k,it) + βU (it) -Π (it) 2 × K k=1 D (k,it) (A (it+1) 1 ) T A (it+1) 1 D (k,it) + βI P -1
(2.137)

U (it+1) = P + αA (it+1) 1 + βA (it+1) 2 + Π (it) 1 + Π (it) 2 /(α + β) (2.138) D (k,it+1) = Diag A (it+1) 1 A (it+1) 2 vec (C (k) ) (2.139) Π (it+1) 1 = Π (it) 1 + γα(A (it+1) 1 -U (it+1) ) (2.140) Chapter 2 2.5. The JDC + ADMM-3 algorithm Π (it+1) 2 = Π (it) 2 + γβ(A (it+1) 2 -U (it+1) ) (2.141)
where I P is a (P × P ) identity matrix, where P + {•} sets all negative values of the input matrix to zero, where Diag(•) returns a diagonal matrix whose diagonal elements are given by the input vector, where denotes the pseudo-inverse operator, where denotes the Hadamard product and denotes the Khatri-Rao product, where vec(•) reshapes a matrix into a column vector by stacking its columns vertically, and where γ ∈ [0, 1.618] [START_REF] Xu | An alternating direction algorithm for matrix completion with nonnegative factors[END_REF]. The pseudo-code of the JDC + ADMM-3 algorithm is summarized in appendix C.7.

Convergence to KKT Points

In this section, following the analysis of the convergence of the ADMM based NMF algorithm [START_REF] Xu | An alternating direction algorithm for matrix completion with nonnegative factors[END_REF], we provide a preliminary study on the convergence of the proposed ADMM algorithm by studying the Karush-Kuhn-Tucker (KKT) conditions of problem (2.129). The KKT conditions of a nonlinear problem are first order necessary conditions for a solution to be optimal, provided that some regularity conditions are satisfied. Consider the following nonlinear optimization problem:

minimize x f (x) (2.142) subject to h i (x) ≤ 0, i = 1, 2, . . . , m g j (x) = 0, j = 1, 2, . . . , r
where x is the optimization variable, where f is the objective function, where h i (i = 1, . . . , m) are the inequality constraint functions, and where g j (j = 1, . . . , r) are the equality constraint functions. Suppose that the objective function f and the constraint functions h i and g j are continuously differentiable at a point x * . We say x * is a KKT point of problem (2.142) if it satisfies the following KKT conditions [START_REF] Boyd | Convex Optimization[END_REF]:

∂f (x * ) ∂x * + m i=1 u i ∂h i (x * ) ∂x * + r j=1 v j ∂g j (x * ) ∂x * = 0 (2.143) h i (x * ) ≤ 0, ∀ i ∈ {1, . . . , m} (2.144)
g j (x * ) = 0, ∀ j ∈ {1, . . . , r} (2.145)
u i ≥ 0, ∀ i ∈ {1, . . . , m} (2.146)
u i h i (x * ) = 0, ∀ i ∈ {1, . . . , m} (2.147)
where u i (i = 1, . . . , m) and v j (j = 1, . . . , r) are constants. Then the KKT point x * is a stationary point with respect to x of the Lagrangian function of problem (2.142).

Regarding the proposed problem (2.129), let us define the following sextuple point:

Z def = (A 1 , A 2 , U , {D (k) }, Π 1 , Π 2 ) (2.148) A point Z is a KKT point of problem (2.

129) if it satisfies the KKT conditions for

Chapter 2 problem (2.129) as follows:

K k=1 (C (k) -A 1 D (k) A T 2 )(A 2 D (k) ) -Π 1 = 0 (2.149) K k=1 (C (k) -A 2 D (k) A T 1 )(A 1 D (k) ) -Π 2 = 0 (2.150) (A 2 A 1 ) T c(k) -(A 2 A 1 )d (k) = 0 (2.151) A 1 -U = 0 (2.152) A 2 -U = 0 (2.153) (Π 1 + Π 2 ) ≤ 0 ≤ U (2.154) (Π 1 + Π 2 ) U = 0 (2.155)
where denotes the Hadamard product and denotes the Khatri-Rao product, where

c (k) = vec(C (k)
), and where d (k) ∈ R P is a column vector containing the diagonal elements of the diagonal matrix D (k) . Then we obtain the following proposition:

Proposition 2.9. Let {Z (it) } ∞ it=1 be a sequence generated by the proposed ADMM algorithm (2.136) -(2.141) that satisfies the following condition:

lim it→∞ (Z (it+1) -Z (it) ) = 0 (2.156)
Then any accumulation point of

{Z (it) } ∞ it=1 is a KKT point of problem (2.129).
The proof of proposition 2.9 is provided in appendix A.7. The following corollary can be derived immediately.

Corollary 2.1. Whenever {Z (it) } ∞ it=1 converges, it converges to a KKT point.

The above preliminary result provides some assurance on the convergence behavior of the ADMM algorithm applied to the non-convex NJDC problem. Further theoretical studies will be considered in future work.

Numerical complexity analysis

This section is devoted to compute the numerical complexities of the proposed methods. The numerical complexity is analyzed in terms of the number of floating point operations (flops). A flop is defined as a multiplication followed by an addition. In practice, only the number of multiplications, required to identify the joint transformation matrix

A ∈ R N ×P + from K symmetric matrices C (k) ∈ R N ×N , k ∈ {1, .
. . , K}, is considered; which does not affect the order of magnitude of the numerical complexity [START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF], Albera et al., 2010, Albera et al., 2012]. The computational complexities of the proposed methods are given in table 2.1 and explained with more details afterwards. In the section of numerical simulations, the proposed methods are compared with five classical JDC methods: the Alternating Columns and Diagonal Center (ACDC) algorithm [Yeredor, 2002], the Fast Frobenius DIAGonalization (FFDIAG) algorithm [Ziehe et al., 2004], the Coupled LU factorization based algorithm (CLU) [Maurandi et al., 2013], and two algorithms based on the general LU factorization and the minimization of two indirect criteria, namely LUJ1D and LUJ2D [Afsari, 2006]. Regarding the CLU method [Maurandi et al., 2013], the JDC criterion (1.27) is approximated by a sequence of (2 × 2) sub-criteria. Therefore, in each Jacobi-like iteration, one elementary lower triangular matrix and one elementary upper triangular matrix can be computed together by EVD, leading to a lower numerical complexity. Furthermore, a NJDC method based on a square change of variable and Levenberg Marquardt optimization [Coloigner et al., 2014c], namely LM + sym , is also included in the comparison. Their complexities are also given in table 2.1.

The JDC +

LU-3 algorithm

The computation of the coefficient matrix G in the cost function (2.36) requires (33N 2 + 7N )K flops. Computing the parameter t i,j by calculating the roots of the 8-th degree polynomial function (2.36) needs 8 3 = 512 flops. Nevertheless, if t i,j is obtained by minimizing the unconstrained cost function (2.44), which is a 4-th degree polynomial function, it costs (15N 2 +4N )K +64 flops. Hence, the complexity of the AC phase, which includes P (N -1) parameters, is between ((15N 2 + 4N )K + 64)P (N -1) and ((33N 2 + 7N )K + 512)P (N -1) flops, depending on the number of unconstrained optimizations is performed. The DC phase consists of estimating K diagonal matrices by means of equation (2.41), and it involves N 2 P K + N 2 P 2 + N 2 P + N P 2 + P 3 + P 2 flops. Hence, the magnitude of the total numerical complexities of one AC phase and one DC phase is O(N 3 P K). The full expression of the complexity of the JDC + LU-3 algorithm is given in table 2.1.

The iJDC +

LU-1 and iJDC + QR-1 algorithms

The iJDC + LU-1 and iJDC + QR-1 algorithms only work on a square matrix A. Therefore, here we assume that N = P . For both algorithms, the inverses of the K matrices C (k) cost N 3 K flops, the initialization of C (k) in equations (2.51) and (2.61) requires 2N 3 K flops. Regarding the estimation of the elementary triangular matrices, in each Jacobilike iteration, the calculation of the parameter t i,j by minimizing the cost function (2.53) needs (5N 2 + 12N -8)K + 64 flops. The calculation cost of A (new) , B (new) and C (k,new) , with k ∈ {1, • • • , K}, is 4N +(4N +1)K flops. Otherwise, if the parameter t i,j is computed without the nonnegativity constraint by means of equation (2.72), it requires 2(N -1)K flops. In this situation, the computation load of A (new) and

C (k,new) is 2N K + N flops.
Regarding the estimation of the Given rotation matrices, in each Jacobi-like iteration, the complexity of calculating the parameters θ i,j by minimizing the cost function (2.65) is equal to (5N 2 +3N +29)K+3375 flops, and the estimation of A (new) , B (new) and C (k,new) , with k ∈ {1, • • • , K}, costs 10N + (12N + 20)K flops. Otherwise, if the parameter θ i,j is obtained by minimizing the unconstrained cost function (2.73), it costs 2K + 8 flops. Then the update of A (new) and C (k,new) needs 4N K + 2N flops. In addition, the row balancing scheme defined in equation (2.74) requires 3N 2 K + N 2 flops. Therefore, the magnitude of global numerical complexities of iJDC + LU-1 and iJDC + QR-1 for estimating all the N (N -1) parameters is between O(N 3 K) and O(N 4 K). Their detailed expressions are given in table 2.1.

Chapter 2

Computational complexity (flops)

(N, P ): the dimension of the joint transformation matrix A; K: the number of the (N × N ) symmetric matrices C (k) ; J i , i ∈ {1, 2, 3, 4, 5, 6, 12}: the number of sweeps (or iterations) used in the i-th method; J rb i , i ∈ {4, 5, 7, 8, 9, 10}: the number of row balancing schemes used in the i-th method; J + i and J * i , i ∈ {7, . . . , 11}: the number of sweeps with and without the nonnegativity constraint used in the i-th method, respectively; CLU, FFDIAG, LUJ1D, LUJ2D, iJDC + LU-1 and iJDC + QR-1 only work on the N = P case.

N = P N > P ACDC (13/3N 3 K + 3N 4 + 2N 2 K + N 3 + N 2 )J 1 (13/3N 2 P K + 3N 3 P + 2N P K + N P 2 + P 2 )J 1 CLU (2N 3 K + 2.5N 2 K + 1.5N 3 + 18.5N 2 -4.5N K -20N )J 2 FFDIAG (2N 3 K + N 3 + 2N 2 K + 4N 2 -4N )J 3 LUJ1D (4N 3 K -6N 2 K + N 3 -N 2 + 2N K)J 4 + (3N 2 K + N 2 )J rb 4 LUJ2D (2N 3 K + 6N 2 K + N 3 + 63N 2 -8N K -64N )J 5 + (3N 2 K + N 2 )J rb 5 LM + sym ( h (N + K) 3 N 3 + h (N + K) 2 N 2 + 4N 6 + N 5 + N 5 K + 11N 4 + 3K 2 N 2 + 3N 3 K + 8N 3 + 3N 2 K + N 2 )J 6
, where h is the number of iterations computing the modification of the Hessian

( h (N + K) 3 P 3 + h (N + K) 2 P 2 + N 4 P 2 + 2N 2 P 4 + N 3 P 3 + N 3 P 2 + N 2 P 3 K + 8N 2 P 2 + 3N P 3 + 3K 2 P 2 + N KP 2 + 2N 2 KP + 6N 2 P + 2N P 2 + 2KP 2 + N KP + P 2 )J 6 iJDC + LU-1 3N 3 K + (4N 3 K -6N 2 K + N 3 -N 2 + 2N K)J * 7 + (5N 4 K + 11N 3 K -23N 2 K + 4N 3 + 60N 2 + 7N K -64N )J + 7 + (3N 2 K + N 2 )J rb 7 iJDC + QR-1 3N 3 K + (4N 3 K -4N 2 K + 1.5N 3 + 2.5N 2 -4N )J * 8 + (5N 4 K + 15.5N 3 K + 5.5N 2 K + 7N 3 + 1712.5N 2 -21N K -1719.5N )J + 8 + (3N 2 K + N 2 )J rb 8 JDC + LU-1 (4N 3 K -6N 2 K + N 3 -N 2 + 2N K)J * 9 + (N 5 + 25N 4 K -22N 3 K -2N 2 K + 4N 3 + 724N 2 + N K -729N )J + 9 + (3N 2 K + N 2 )J rb 9 (4N 2 P K -6N P K + N 2 P -N P + 2P K)J * 9 + (N 4 P + 25N 3 P K + 2N 3 K - 24N 2 P K -2N P K + 4N 2 P + 724N P + P K -729P )J + 9 + (3N 2 K +N P )J rb 9 JDC + LU-2 (2N 3 K + 6N 2 K + N 3 + 63N 2 -8N K -64N )J * 10 + (5N 5 K + 29N 4 K + 2N 5 - 16N 3 K + 3N 4 -19N 2 K + 3N 3 + 1323N 2 + 3N K -1331N )J + 10 + (3N 2 K + N 2 )J rb 10 (2N 2 P K + 6N P K + N 2 P + 63N P -8P K -64P )J * 10 + (5N 4 P K + 29N 3 P K + 2N 4 P + 2N 3 K -18N 2 P K + 3N 3 P -19N P K + 3N 2 P + 1323N P + 3P K - 1331P )J + 10 + (3N 2 K + N P )J rb 10 JDC + LU-3 (15N 4 K -10N 3 K + N 4 -4N 2 K + 3N 3 + 65N 2 -64N )J * 11 + (33N 4 K - 25N 3 K + N 4 -7N 2 K + 3N 3 + 513N 2 -512N )J + 11 (15N 3 P K -10N 2 P K + N 2 P 2 -4N P K + N 2 P + N P 2 + P 3 + P 2 + 64N P - 64P )J * 11 + (33N 3 P K -25N 2 P K + N 2 P 2 -7N P K + N 2 P + N P 2 + P 3 + 512N P + P 2 -512P )J + 11 JDC + ADMM-3 (3N 3 K + N 4 + 6N 2 K + 7N 3 + 8N 2 + 2N )J 12 (3N 2 P K + N 2 P 2 + 2N P K + 4P 2 K + N 2 P + 3N P 2 + 3P 3 + 7N P + P 2 + 2P )J
Table 2.1: Numerical complexities of the proposed methods and of ACDC [Yeredor, 2002], CLU [Maurandi et al., 2013], FFDIAG [Ziehe et al., 2004], LUJ1D, LUJ2D [Afsari, 2006] and LM + sym [Coloigner et al., 2014c].

Chapter 3

Numerical simulations

In this chapter, the performance of the proposed algorithms is evaluated through simulated data fitting the NJDC model. Several experiments are designed to study the convergence property, the influence of SNR, the impact of the number K of the input symmetric matrices C (k) , the effect of the coherence of the column vectors of the joint mixing matrix A and the influence of the modulus of uniqueness of JDC. In each experiment, the performance is analyzed for both estimating a square joint transformation matrix A and a non-square matrix A. The proposed algorithms are compared with five classical nonorthogonal JDC methods, which do not exploit the nonnegativity property of A, namely the Alternating Columns and Diagonal Center (ACDC) algorithm [Yeredor, 2002], the Fast Frobenius DIAGonalization (FFDIAG) algorithm [Ziehe et al., 2004], the Coupled LU factorization based algorithm (CLU) [Maurandi et al., 2013], and two algorithms based on the general LU factorization and the minimization of two indirect criteria, namely LUJ1D and LUJ2D [Afsari, 2006]. In addition, a NJDC method based on a square change of variable and Levenberg Marquardt optimization [Coloigner et al., 2014c], namely LM + sym , is also included in the comparison.

Experimental setup

K synthetic symmetric matrices C (k) ∈ R N ×N with k ∈ {1, . . . , K} are generated randomly according to the NJDC model (1.51). When used without further specification, all the algorithms are manipulated under the following conditions:

1. Model generation: the joint transformation matrix A ∈ R N ×P + is randomly drawn from a uniform distribution on the interval [0, 1]. The diagonal elements of the K diagonal matrices D (k) with k ∈ {1, . . . , K} are drawn from the zero-mean unitvariance Gaussian distribution. Each pure matrix C (k) is perturbed by a random symmetric noise matrix as follows:

C (k) N = C (k) C (k) F + σ (k) R (k) (R (k) ) T R (k) (R (k) ) T F (3.1)
where σ (k) is a scalar controlling the noise level, and where the components of the noise matrix R (k) ∈ R N ×N are drawn from the zero-mean unit-variance Gaussian distribution. Then the SNR is defined by SNR = -20 log 10 (σ (k) ).

Convergence property

Chapter 3

2. Initialization: in each Monte Carlo trial, all the algorithms are initialized by the absolute value of the inverse of the output of the Uniformly Weighted Exhaustive Diagonalization with Gauss itErations (UWEDGE) algorithm [START_REF] Tichavský | Fast approximate joint diagonalization incorporating weight matrices[END_REF] due to its fast speed.

3. Afsari's row balancing scheme: the LUJ1D, LUJ2D, iJDC

+ LU-1 , iJDC + QR-1 , JDC + LU-1
and JDC + LU-2 algorithms perform the row balancing scheme once per each run of five sweeps.

4. Stopping criterion: all the algorithms stop either when the relative error of the corresponding criterion between two successive sweeps of a Jacobi-like method (or iterations of an alternating minimization method) is less than 10 -4 or when the number of sweeps (iterations) exceeds 500.

5. Performance measurement: the performance is measured by means of the error between the true joint transformation matrix A and the estimate A, and the numerical complexity. Note that the numerical complexity can be approximately measured by computing the CPU time. However, the CPU time is computer-dependent and also varies according to the operating system. In this report, the simulations are carried out in Matlab v7.14 on Mac OS X and run on Intel Quad-Core CPU 2.8 GHz with 32 GigaByte (GB) memory.

In order to assess the estimation accuracy of A, we use the following scale-invariant and permutation-invariant distance [Coloigner et al., 2014a,Coloigner et al., 2014c]:

α(A, A) = 1 N N n=1 min (n,n )∈I 2 n d(a n , a n ) (3.2)
where a n and a n are the n-th column vector of A and the n -th column vector of A, respectively. I 2 n is a set of indices defined recursively by

I 2 1 = {1, • • • , N } × {1, • • • , N }, and I 2 n+1 = I 2 n -J 2 n where J 2 n = argmin (n,n )∈I 2 n d(a n , a n ).
In addition, d(a n , a n ) is the pseudo-distance between two column vectors defined by [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF]:

d(a n , a n ) = 1 - a T n a n 2 a n 2 a n 2 (3.3)
The criterion (3.2) is an upper bound of the optimal permutation-invariant criterion. It avoids the burdensome computation of all the permutations. A small value of (3.2) means a good performance in the sense that A is close to A. Moreover, we repeat all the experiments with 500 Monte Carlo trials.

Convergence property

In this section, we aim at studying the convergences of the proposed algorithms, as well as those of the ACDC, CLU, FFDIAG, LUJ1D, LUJ2D and LM + sym algorithms.

Test on square matrices

In this experiment, the dimensions of the square joint transformation matrix A ∈ R N ×P + are set to N = 3 and P = 3. The number of the input matrices C (k) is fixed to K = 5.

Chapter 3

Convergence property

The convergence is assessed under three SNR conditions: SNR = 30 dB, 10 dB and 0 dB, respectively. Figure 3.1 shows the average convergence curves measured in terms of the cost functions as a function of sweeps/iterations. The twelve algorithms are separated into four groups according to their cost functions. In general, each algorithm exhibits a relative stable convergence pattern under different SNR levels. Among the classical approaches, CLU shows the fastest convergence property. Generally, it converges in less than 10 sweeps. FFDIAG also converges fast when SNR is equal to 10 dB and 0 dB. LUJ1D and LUJ2D require less than 20 sweeps to converge. The alternating optimization based method ACDC, as well as LM + sym , converge slower than those Jacobi-like algorithms. Generally ACDC converges after 200 iterations, while LM + sym requires more than 500 iterations to converge. Regarding the proposed methods, JDC + LU-1 converges in 30 sweeps. JDC + LU-2 converges in 20 sweeps. iJDC + LU-1 and iJDC + QR-1 converge in less than 40 sweeps and 100 sweeps, respectively. Generally, JDC + ADMM-3 converges in 50 iterations. JDC + LU-3 converges the slowest among the proposed methods and requires more than 500 iterations to converge. In addition, it can be seen that there exist some fluctuations in the convergence curves of CLU, FFDIAG, JDC + LU-2 and JDC + ADMM-3 , which are caused by the increase of the cost functions in a few Mont Carlo realizations. For CLU and FFDIAG, they minimize the approximations of the cost function (1.27). Therefore, in a few cases, the inaccurate estimation of the parameters leads to an increase of the original cost function. Similarly, the fluctuation of the convergence curve of JDC + LU-2 is caused by the step of unconstrained minimization (2.127), which is also an approximation of the original JDC criterion. Regarding JDC + ADMM-3 , the limited fluctuation is due to the projection of U into the nonnegative space by equation (2.138), which can lead to an increase of the cost function in rare cases [START_REF] Kim | Fast Newton-type methods for the least squares nonnegative matrix approximation problem[END_REF]. In most cases, these algorithms show good convergence properties. The numerical complexities of all the methods in this experiment are listed in table 3.1. It shows that the classical Jacobi-like methods, CLU, FFDIAG, LUJ1D and LUJ2D achieve the lowest numerical complexities. The proposed methods require more flops to converge, but they are still much less costly than LM + sym . Among the proposed methods, JDC + ADMM-3 is the cheapest method in terms of complexity, and it needs just a bit more flops than ACDC.

Test on non-square matrices

In this experiment, we compare four of the proposed algorithms, namely JDC + LU-1 , JDC + LU-2 , JDC + LU-3 and JDC + ADMM-3 , with ACDC and LM + sym for estimating a non-square matrix A. These methods can be directly applied to solve a non-square JDC problem.

The dimensions of the joint transformation matrix A ∈ R N ×P + are set to N = 5 and P = 3, and the number of the input matrices C (k) is fixed to K = 5. The convergence is also assessed under three SNR conditions: SNR = 30 dB, 10 dB and 0 dB, respectively. The convergence curves of those six algorithms measured in terms of the cost functions as a function of sweeps/iterations are shown in figure 3.2. It shows that JDC + LU-1 and JDC + LU-2 converge in less than 50 sweeps and 10 sweeps, respectively, whatever the SNR is. ACDC and JDC + ADMM-3 show similar convergence behaviors, and they require 50 iterations to converge. LM + sym and JDC + ADMM-3 converge in 200 iterations when SNR is equal to 0 dB, and in 500 iterations when SNR is equal to 10 dB. Nevertheless, it seems that they do not converge in 500 iterations when SNR = 30 dB. The numerical complexities of those six methods are listed in table 3.2. We can observe that ACDC and JDC + ADMM-3 maintain the lowest numerical complexities. LM + sym requires the heaviest computational load when SNR is equal to 10 dB and 0 dB. When SNR = 30 dB, JDC + LU-3 is the most computationally expensive method. The two Jacobi-like methods, JDC + LU-1 and JDC + LU-2 need relatively moderate complexities, while the latter is more efficient than the former for estimating a non-square matrix.

Influence of SNR

In this section, we study the behaviors of the considered algorithms as a function of SNR.

Test on square matrices

In this experiment, the dimensions of the square joint transformation matrix A ∈ R N ×P + are fixed to N = 3 and P = 3. The number of the input matrices C (k) is set to K = 5. We repeat the experiment with SNR ranging from -20 dB to 40 dB with a step of 2 dB. Figure 3.3(a) depicts the average curves of the estimation error α(A, A) of the twelve algorithms as a function of SNR. The obtained results show that the performance of all the methods increases quasi-quadratically as SNR grows. For the classical unconstrained methods, ACDC performs better than CLU, FFDIAG, LUJ1D and LUJ2D, while the latter four Jacobi-like methods show almost the same performance. The nonnegativity constraint obviously helps LM + sym and the proposed six methods to improve the estimation accuracy especially for lower SNR values. When SNR is less than 20 dB, the proposed four Jacobi-like algorithms, iJDC + LU-1 , iJDC + QR-1 , JDC + LU-1 and JDC + LU-2 outperform the others, while the two alternating minimization based NJDC methods JDC + LU-3 and JDC + ADMM-3 , as well as LM + sym , also product better results than the classical unconstrained algorithms. Figure 3.3(b) displays the detail information for some selected methods with SNR ranging from 20 dB to 40 dB. We can seen that JDC + ADMM-3 maintains a competitive advantage over the other methods under those higher SNR conditions. JDC + LU-3 also achieves good results, and it outperforms ACDC when SNR ≤ 28 dB. LM + sym , JDC + LU-1 and JDC + LU-2 give similar performance, and they produce more accurate estimation results than iJDC + LU-1 and iJDC + QR-1 . The average numerical complexities at each SNR level of all the methods in this experiment are shown in figure 3.3(c). It shows that the classical Jacobi-like methods CLU, FFDIAG and LUJ1D obtain the lowest numerical complexities. Generally, LM + sym is the most computationally expensive algorithm. Among the proposed methods, the complexity of iJDC + LU-1 is the lowest, and increase as SNR grows. The complexity of JDC + ADMM-3 is almost identical to that of ACDC. iJDC + QR-1 , JDC + LU-1 , JDC + LU-2 and JDC + LU-3 cost more flops to achieve better estimations of A, but they are still much more efficient than LM + sym . Therefore, the JDC + LU-2 and JDC + ADMM-3 give the best performance/complexity trade-off for all the considered SNR values.

Test on non-square matrices

In this experiment, a non-square joint transformation matrix A ∈ R N ×P + is considered, where the dimensions are set to N = 5 and P = 3, and where the number of the input matrices C (k) is fixed to K = 5. We repeat the experiment with SNR ranging from -20 dB to 40 dB with a step of 2 dB. Figure 3.4(a) shows the average curves of the estimation error α(A, A) as a function of SNR for ACDC, LM + sym , JDC + LU-1 , JDC + LU-2 , JDC + LU-3 and JDC + ADMM-3 . As illustrated in the picture, the performance of all these methods increases as SNR grows. The nonnegativity constrained algorithms, LM + sym , JDC + LU-1 , JDC + LU-2 , JDC + LU-3 and JDC + ADMM-3 , outperform ACDC when SNR is less than 14 dB, and perform similarly to ACDC when SNR exceeds 20 dB. The proposed JDC + LU-2 algorithm is more superior than the other methods when SNR is less than 12 dB. More detailed information of all the methods with SNR ranging from 12 dB to 40 dB is displayed in figure 3.4(b). Under these higher SNR levels, JDC + ADMM-3 gives the most accurate estimation results. JDC + LU-3 outperforms ACDC when SNR < 26 dB. JDC + LU-1 and JDC + LU-2 also achieve comparable results. The average numerical complexity over all the SNR values for each of the estimation error α(A, A) of the ACDC, LM + sym , JDC + LU-1 , JDC + LU-2 , JDC + LU-3 and JDC + ADMM-3 algorithms as a function of K are displayed in figure 3.6(a). The result shows that the increase of K benefits most JDC algorithms. JDC + LU-2 gives the most accurate estimation results when only a small number of input matrices is available (e.g. K ≤ 7). When we have a large number of input matrices, that is to say K is larger than 11, JDC + ADMM-3 maintains competitive advantages over the other methods with large performance gains. LM + sym , JDC + LU-2 and JDC + LU-3 outperform ACDC when K is less than 15, and show similar performance to ACDC when K > 15. The result of JDC + LU-1 is less desirable for estimating a non-square matrix A when K is large. Figure 3.5(b) depicts the average numerical complexities of all the algorithms in this experiment. The complexities of all the methods in terms of flops grows as K increases. Obviously, LM + sym requires the highest numerical complexity, while ACDC and JDC + ADMM-3 are the most computationally efficient methods. The complexity of JDC + LU-2 is nearly one order of magnitude higher than that of ACDC. The complexities of JDC + LU-1 and JDC + LU-3 are between those of JDC + LU-2 and LM + sym . Consequently, the JDC + ADMM-3 algorithm provides the best performance/complexity compromise, while the JDC + LU-2 and JDC + LU-3 algorithms are also competitive with regard to LM + sym .

Influence of coherence of matrix A

In this section, the impact of the coherence between two column vectors of the joint transformation matrix A on the performance of the considered algorithms is evaluated. A nonnegative matrix A of dimension (N × P ) is generated with a coherence parameter β as follows:

A = [a 1 , (1 -β)a 1 + βa 2 ,

Influence of modulus of uniqueness

In this section, the effect of the modulus of uniqueness ρ of a JDC problem is evaluated. The modulus of uniqueness ρ is defined in equation (1.26) in section 1.3.2. By its definition, ρ falls in the range of [0, 1]. The JDC problem is considered to be ill-conditioned when ρ is close to 1. In order to control ρ, firstly we randomly generate an orthogonal matrix D ∈ R K×P so that ρ = 0 by orthogonalizing a (K × P ) random matrix whose components are drawn from the zero-mean unit-variance Gaussian distribution. Secondly we rotate the P column vectors of D such that all the internal angles between any column vectors are equal to a predefined value ψ. Therefore ρ is only controlled by the angle ψ and is equal to | cos(ψ)|. Finally we can obtain K diagonal matrix D (k) ∈ R P ×P , where the diagonal elements of each D (k) are equal to the elements of the k-th row vector of the predefined matrix D.

Test on square matrices

In this simulation, the dimensions of the square joint transformation matrix A ∈ R N ×P + are set to N = 3 and P = 3, the number of the input matrices C (k) is fixed to K = 5, and the SNR value is set to 10 dB. We repeat the experiment with the angle ψ ranging from 0 to π/2 with a step of π/60. A small ψ value means a large ρ value. Figure 3.9(a) shows the average curves of the estimation error α(A, A) of the twelve algorithms as a function of ψ. As we can see, in general cases, the performance of all the algorithms increases as the value of ψ grows. Four classical Jacobi-like methods CLU, FFDIAG, LUJ1D and LUJ2D give moderate results, while ACDC offers an improvement in terms of estimation accuracy. The nonnegativity constraint obviously benefits LM + sym and the proposed six methods to further enhance the results. Generally speaking, the performance of the proposed four Jacobi-like NJDC methods, namely JDC + LU-1 , JDC + LU-2 , iJDC + LU-1 and iJDC + QR-1 , surpasses that of LM + sym , JDC + LU-3 and JDC + ADMM-3 , where the latter three

29). The JDC +

LU-1 algorithm maintains competitive advantage through all the considered ψ values. The average numerical complexities of all the algorithms as a function of ψ in this experiment is shown in figure 3.9(b). As illustrated in the figure, generally the classical Jacobi-like methods CLU and FFDIAG preserve the lowest numerical complexities. LUJ1D and LUJ2D are still efficient. The complexity of iJDC + LU-1 is the lowest among those of the proposed six methods, and it is also lower than or equal to that of ACDC. JDC + ADMM-3 costs a bit more computation load than ACDC, followed by JDC + LU-2 and iJDC + QR-1 , consecutively. The complexities of JDC + LU-1 and JDC + LU-3 are almost one order of magnitude higher than that of JDC + ADMM-3 . The LM + sym algorithm requires the heaviest workload. Consequently, the proposed iJDC + LU-1 and JDC + LU-2 algorithms offer the best performance/complexity trade-off for all the considered ψ values. value of modulus of uniqueness is involved, however, at the cost of a higher computational complexity. When a well-conditioned NJDC problem is considered, that is to say, when a high value of SNR can be expected or a large number of input matrices is available, the JDC + ADMM-3 algorithm is the most effective in terms of estimation accuracy and also has a low numerical complexity.

Chapter 4

Blind source separation applications

We have shown the good performance of the proposed algorithms for solving simulated NJDC model in the last chapter. The purpose of this chapter is to illustrate the practical usefulness of our methods by solving three BSS problems from very different application domains. The first one which is rooted in the biomedical diagnostics, consists of separating synthetic realistic Magnetic Resonance Spectroscopy (MRS) signal. The second one is an automatic music transcription problem. It aims at estimating the musical notes and their attributes from a single-channel music record by decomposing its magnitude spectrogram. The third one is devoted to separating digital images which are degraded by the so-called show-through effect. Such a problem is usually caused by the seeping of ink from the reverse side or scanning a semi-transparent paper. In the following paragraphs, the SeNICA methods based on the proposed NJDC algorithms, namely, for example, JDC + ADMM-3 -ICA, consist of the following four steps:

1. Compressing the M realizations of the observation vector x ∈ R N by means of a nonnegative matrix W ∈ R P ×N + computed using the NonNegative COMPression (NN-COMP) method, which is addressed in appendix B, leading to a compressed observation vector x ∈ R P , where N , P and M are the numbers of the observations, of the sources and of the samples points, respectively. The numerical complexity of NN-COMP is analyzed in equation (B.20) of appendix B.

Estimating the FO cumulant array C

(4)

x ∈ R N ×N ×N ×N from the original observation x, or estimating the FO cumulant array C (4)

x ∈ R P ×P ×P ×P from the compressed observation x; then extracting all the cumulant matrix slices to form a matrices set. The complexity of computing C (4)

x is given as follows [START_REF] Albera | Biomedical applications[END_REF]:

f (4) cum (N, M ) = 3 8 N (N + 1)(N 2 + N + 2)M (4.1)
while the complexity of calculating C

(4)

x is f x , respectively. The numerical complexities of the proposed NJDC methods are listed in table 2.1.

Chapter 4

give good results, but their separated spectra are still disturbed. Regarding NMF and semi-NMF, the Choline peak still slightly exists in the estimated Myo-inositol spectrum, and the estimated Choline spectrum is also disturbed. SOBI cannot separate both source spectra. These experimental results show that by fully exploiting the prior information of data, such as the nonnegativity of the mixing matrix and the statistical independence of the sources, the proposed methods achieve visually better estimation results. In the following sections, the performances of all the considered methods are quantitatively analyzed by means of two experiments. More precisely, we study the influences of the number of observations and of the SNR.

Influence of the number of observations

In this experiment, the effect of the number of observations N is evaluated. The SNR is fixed to 10 dB. The fifteen methods are compared with N ranging from 3 to 41. The average curves of the estimation error α(s T , s T ) of all the methods as a function of N are shown in figure 4.4(a). It can be seen that the estimation errors of all the methods improves as N increases. It suggests that in noisy BSS contexts, using more sensors often yields better results [START_REF] Joho | Overdetermined blind source separation: Using more sensors than source signals in a noisy mixture[END_REF]. Obviously, the proposed ten methods maintain the competitive advantages when N is greater than 5, where the six methods using the NN-COMP preprocessing slightly outperform the ones without a compression step. We can infer that NN-COMP not only compresses the dimension but also reduces the effect of additive noise to a certain extent. Moreover, in a relatively low SNR context, the merit of the NN-COMP step overweights its disadvantage, for instance, introducing additional numerical errors as in the classical prewhitening step. Regarding the classical BSS methods, CoM2 and NICA behave quite similarly and provide good results when N is less than 25. However, they restrict the search space to the group of orthogonal matrices, which may limit their BSS performance in this problem. As we can see, their performance stabilizes when N exceeds 20. NMF outperforms CoM2 and NICA when N ≥ 30. Semi-NMF gives less accurate estimation results than NMF. It seems that only the nonnegativity constraint cannot guarantee a good separation result. SOBI yields the worst results. It may be because the source processes in this experiment do not fully satisfy the non-iid and stationary assumption of SOBI. The average curves of the numerical complexities of all the considered methods in this experiment are shown in figure 4.4(b). We can notice that the numerical complexities of all the methods increase as a function of N . CoM2 is the most computationally efficient method. The proposed six methods incorporating NN-COMP cost a bit more flops than CoM2, where NN-COMP + JDC + LU-3 -ICA is somewhat more expensive than the other five algorithms. NICA and SOBI are still economical in terms of flops. The aforementioned nine methods keep relatively low numerical complexities since all of them use a dimension reduction step by means of either classical prewhitening or NN-COMP. The methods without using a dimension reduction step, namely NMF, semi-NMF, JDC + LU-1 -ICA, JDC + LU-2 -ICA, JDC + LU-3 -ICA and JDC + ADMM-3 -ICA, are more computationally intensive. NMF requires the heaviest workload when N is less than 20. The complexity of semi-NMF is almost half of that of NMF. The proposed four methods JDC + LU-i -ICA, i ∈ {1, 2, 3} and JDC + ADMM-3 -ICA preserve moderate complexities when a small value of N is considered, say N < 15. However, their complexities show a higher growth rate as N increases, and surpass that of NMF when N is above 22. The large amount of computation load of these four methods partially comes from the estimation of the FO cumulant arrays. Among these four algorithms, JDC + LU-3 -ICA seems the most costly, while JDC + ADMM-3 -ICA is less computationally expensive. Therefore, the proposed SeNICA methods incorporating the NN-COMP offer the best performance/complexity compromise in this experimental context.

Influence of SNR

In this simulation, we study the influence of SNR on the performance of all the considered fifteen methods. The number of observations N is set to 20. SNR is varied from 0 dB to 60 dB with a step of 2 dB. The average curves of the estimation error α(s T , s T ) of all methods as function of SNR are presented in figure 4.5(a). It show that the performance of all methods increases as SNR grows, and then practically stabilizes for high values of SNR. Obviously, the proposed SeNICA algorithms offer the most accurate separation solutions, especially when SNR is above 5 dB. The performance of the proposed ten methods is Chapter 4 almost identical when the value of SNR is between 0 dB and 18 dB, where the six methods using NN-COMP are slightly more preferable when SNR < 12 dB. As far as higher values of SNR are considered (e.g. SNR > 20 dB), the proposed four methods, namely JDC + LU-i -ICA, i ∈ {1, 2, 3} and JDC + ADMM-3 -ICA, perform better than the six methods incorporating NN-COMP. This phenomenon suggests that in a high SNR context, the additional numerical errors introduced by NN-COMP may bring some limitations for the ensuing SeNICA algorithms. Regarding the classical BSS methods, when SNR is lower than 20 dB, CoM2 and NICA offer good results, followed by NMF, semi-NMF and SOBI, consecutively. The performance of CoM2, NICA, NMF and semi-NMF is almost the same when SNR is above 24 dB. SOBI achieves the worst solutions. Figure 4.5(b) illustrates the average curves of the numerical complexities of all methods in this experiment. We can observed that the complexity evolution is quite stable with the change in SNR. CoM2 maintains the lowest complexity, followed by the proposed six methods combined with NN-COMP, where NN-COMP + JDC + LU-3 -ICA is more expensive. NICA and SOBI are also efficient. NMF is the most computationally intensive. Semi-NMF costs nearly half workload than NMF. The complexities of the proposed four methods without using NN-COMP are between those of semi-NMF and NMF. Consequently, the proposed SeNICA methods incorporating NN-COMP provide the best performance/complexity tradeoff in this experiment. Without using NN-COMP, the proposed four methods JDC + LU-i -ICA, i ∈ {1, 2, 3} and JDC + ADMM-3 -ICA give the most accurate estimation results at a cost of higher but acceptable numerical complexities, when a high SNR value is provided.

The results of the above two experiments demonstrate the improvement of the proposed methods in terms of the source estimation accuracy, and also show that exploiting the two a priori of the data, namely the nonnegativity of the mixing matrix and the statistical independence of the sources, allows us to achieve better estimation results.

Automatic music transcription

Estimating the musical notes and their attributes (pitch, duration and onset time) from a music record is an objective of automatic music transcription. To this end, many unsupervised algorithms propose to decompose the music magnitude spectrogram X as X = AS. Such a linear decomposition consists in modeling a time-frequency representation of music signal as a sum of basic atoms S. Each row vector of matrix S represents a basic atom, which corresponds to the harmonically fixed spectrum of a single note [Smaragdis andBrown, 2003, Fuentes et al., 2013]. An atom will be active whenever the corresponding note is played. Furthermore, each column vector of matrix A is modeled as the nonnegative temporal weights of a basic note spectrum, that describes the onset time and duration of an active note. Many algorithms based on ICA [START_REF] Brown | Independent component analysis for onset detection in piano trills[END_REF], Plumbley and Abdallah, 2003, Brown and Smaragdis, 2004], NICA [Plumbley, 2003, Dittmar andUhle, 2004] and NMF [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF], Holzapfel and Stylianou, 2008, Févotte et al., 2009, Gao et al., 2011, Smaragdis et al., 2014] were proposed to solve such a problem. Nevertheless, the temporal weights containing negative values obtained by classical ICA methods bring some difficulties for the physical interpretation.

In this experiment, a single-channel polyphonic piano record played from the music sheet in figure 4.6(a) is adopted [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: with application to music analysis[END_REF]. Such a 15.6 seconds record sampled at 22.05 kHz is made up from four different notes: C 5 , A b 4 , F 4 and D b 4 (from top to bottom in figure 4.6(a)). Theoretically a note consists of several related sine waves called the fundamental and the harmonics. The sound of a note has a harmonic frequency spectrum [Burg, 2008]. The lowest frequency present is the fundamental frequency, and is the frequency at which the entire wave vibrates. A harmonic is a component frequency of a signal that is an integer multiple of the fundamental frequency. For example, if the fundamental frequency is f , the harmonics have frequencies 2f , 3f , 4f , ... , etc. The theoretical fundamental frequencies of four notes in this experiment are given in table 4.1. The piano record is contaminated by an additive pink noise with a SNR value of 5 dB. A pink noise occurs in many digital recording devices, whose power spectral density is inversely proportional to the frequency of the signal. The spectrogram of the noisy record is computed using a hamming window of 125 milliseconds length with 50% overlap between two frames, leading to 251 frames and 1379 frequency bins. Moreover, the time domain waveform of the record and its magnitude spectrogram are presented in figures

Document restoration

Some digital images of documents are degraded by the so-called show-through effect, which is usually caused by the seeping of ink from the reverse side or scanning a semitransparent paper [START_REF] Tonazzini | Multichannel blind separation and deconvolution of images for document analysis[END_REF]. In [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF], Tonazzini et al., 2007], Tonazzini et al. modeled the show-through effect as a linear instantaneous superimposition of the back and the front sides of the scanned images, and achieved the separation of two sides of images by using ICA. Thereafter, Merrikh-Bayat et al. proposed a NMF based technique for removing show-through [START_REF] Merrikh-Bayat | Using non-negative matrix factorization for removing show-through[END_REF]. It suggests that imposing the nonnegativity constraint can improve the result. In addition, better separation result can be achieved by considering more sophisticated models, such as the linear convolutive mixing model [START_REF] Tonazzini | Multichannel blind separation and deconvolution of images for document analysis[END_REF], the linear-quadratic mixing model [START_REF] Merrikh-Bayat | Linear-quadratic blind source separating structure for removing show-through in scanned documents[END_REF] and the nonlinear mixing model [START_REF] Almeida | Nonlinear separation of show-through image mixtures using a physical model trained with ICA[END_REF], however, leading to higher numerical complexities.

In this experiment, we apply six BSS methods to a real double-sided grayscale document of dimension (1136 × 1038), which is suffered from strong show-through effect. The considered methods include CoM2 [Comon, 1994], NICA [Plumbley, 2003], a projected gradient based NMF [START_REF] Merrikh-Bayat | Using non-negative matrix factorization for removing show-through[END_REF], semi-NMF [Ding et al., 2010], and two of the proposed methods, namely JDC + LU-2 -ICA and JDC + ADMM-3 -ICA. The input front side and back side images displayed in figures 4.9(a) serve as 2 observations. Figures 4.9(b) to 4.9(g) show the restored front side and back side images by the considered six methods, respectively. It can be seen that for CoM2, NMF and the proposed methods, the show-through effect is mostly canceled. Their separation results are similar, although some small discrepancies exist in the detailed structure of the restored images, such as the straight line in the bottom of the left image. The results of JDC + LU-2 -ICA and JDC + ADMM-3 -ICA are visually identical. NICA and semi-NMF seem not able to separate both sides of images. Furthermore, the numerical complexities of all the methods are listed in table 4.3. It shows that NMF requires the heaviest workload, while the proposed JDC + LU-2 -ICA and JDC + ADMM-3 -ICA algorithms are the most computationally efficient. Hence, the proposed methods offer the best performance/complexity compromise in this experiment.

Concluding remark

SeNICA turns up in many real-life BSS problems, which involve nonnegative static mixtures, such as in the blind separation of MRS signal, automatic music transcription of single-channel piano record and restoration of documents degraded by the show-through effect. The proposed SeNICA methods are applied to solve these problems. We compare them to five state-of-the-art BSS algorithms, namely two efficient ICA methods CoM2 [Comon, 1994] and SOBI [Belouchrani et al., 1997], NICA [Plumbley, 2003], NMF [Kim andPark, 2008, Merrikh-Bayat et al., 2010] and semi-NMF [Ding et al., Chapter 4 2010]. Based on the presented studies, NICA, SOBI and semi-NMF sometimes fail to separate the sources. CoM2 and NMF offer better results, despite some imperfections. The proposed methods are most likely to provide better results than the classical methods, in terms of estimation accuracy of the sources, meanwhile maintaining moderate numerical complexities. These BSS applications demonstrate the validity and potential usefulness of the proposed algorithms. The experimental results have also shown that by fully exploiting the prior information of data, such as the nonnegativity of the mixing matrix and the statistical independence of the sources, the proposed methods achieve better estimation results.

Chapter 5

Conclusion and perspectives

Conclusion

The purpose of this thesis is motivated by the SeNICA problem, where the mixing matrix is constrained to have nonnegative components, and where the statistically independent sources may have mixed signs. Such a problem can be solved by means of JDC algorithms with a nonnegativity constraint on the joint transformation matrix, leading to the concept of NJDC. Therefore, we develop six novel NJDC algorithms, which are presented thoroughly in chapter 2. The six proposed methods can be separated into two groups according to the way of performing the optimization. The first group, containing five algorithms, is based on the Jacobi-like optimization. In this group, the nonnegativity constraint is ensured by a square change of variable. The general idea of our Jacobi-like optimization is to factorize the square root of the nonnegative joint transformation matrix as a product of elementary matrices which is defined by only one parameter, then to estimate these elementary matrices one by one in a specific order. The first method adopts the LU matrix factorization and performs the optimization on a simple directfit criterion, however, suffers from a high numerical complexity. In order to reduce the numerical complexity, we reformulate the indirect-fit criterion and consider both the LU and QR factorizations. The proposed reformulation simplifies the optimization and leads to two computationally efficient algorithms. Nevertheless, the validity of these two algorithms relies on some crucial assumptions (A is square and D (k) does not contain zero diagonal elements). Consequently, the fourth and fifth algorithms that depend on less assumption are proposed. These two methods are based on a direct optimization of two indirect-fit criteria incorporating the LU matrix factorization. The second group, containing one algorithm, is based on the ADMM optimization. This algorithm is motivated by the fact that the NJDC problem can be recast as a direct least-squares fit problem subject to symmetry and nonnegativity constraints. Then such an algorithm is derived by successively minimizing the augmented Lagrangian function of the direct-fit criterion with respect to the variables and the multipliers. In addition, the numerical complexity of each algorithm has been analyzed, allowing a fair comparison of their performance.

In chapter 3, the performance of the proposed NJDC algorithms is evaluated on simulated jointly diagonalizable matrices. Five classical nonorthogonal JDC methods without nonnegativity constraints, including ACDC [Yeredor, 2002], CLU [Maurandi 5.1. Conclusion Chapter 5 et al., 2013], FFDIAG [Ziehe et al., 2004], LUJ1D [Afsari, 2006], LUJ2D [Afsari, 2006] and one NJDC method LM + sym [Coloigner et al., 2014c], are tested as baseline methods. The performance is assessed in terms of the matrix estimation accuracy and the numerical complexity. Throughout the chapter, i) the convergence property, ii) the influence of model errors, iii) the impact of the number of input matrices, iv) the effect of coherence of the column vectors of the joint transformation matrix, and v) the influence of the modulus of uniqueness, are extensively studied through Monte Carlo experiments. Simulation results show that a better estimation accuracy can be achieved particularly in difficult contexts by means of exploiting the nonnegativity a priori. For a difficult NJDC problem with either a lower value of SNR, a small number of input matrices, a high level of coherence presented in the joint transformation matrix or a large value of modulus of uniqueness, generally the proposed Jacobi-like algorithms provide the most accurate estimation results. Among them, the JDC + LU-2 algorithm offers the best compromise between estimation accuracy and numerical complexity. When the considered NJDC problem is well-conditioned, that is to say, when a high value of SNR can be expected or a large number of input matrices is available, the algorithm based on the ADMM optimization, namely JDC + ADMM-3 , is the most effective in terms of the estimation accuracy and is also computationally efficient. Compared to an existing nonnegative method based on the Levenberg Marquardt optimization [Coloigner et al., 2014c] whose performance is also competitive in some experiments, the proposed methods require less computation load, and save the total numerical complexity by up to three orders of magnitude. Generally speaking, the proposed JDC + LU-2 and JDC + ADMM-3 algorithms offer the best performance/complexity tradeoff.

In chapter 4, we apply the proposed NJDC algorithms on a set of FO cumulant matrix slices, giving birth to a class of SeNICA methods. SeNICA turns up in quite a number of real-life BSS problems, which involve nonnegative static mixtures, such as in the blind separation of MRS signals, automatic music transcription of single-channel piano records and restoration of documents degraded by the show-through effect. The proposed SeNICA methods are applied to solve these problems in order to verify their source separation ability. Our methods are compared with five state-of-the-art BSS algorithms, namely two efficient ICA methods CoM2 [Comon, 1994] and SOBI [Belouchrani et al., 1997], three nonnegativity constrained methods NICA [Plumbley, 2003], NMF [Kim andPark, 2008, Merrikh-Bayat et al., 2010] and semi-NMF [Ding et al., 2010]. Based on the obtained results, the classical methods, especially NICA, SOBI and semi-NMF sometimes fail to separate the latent sources. CoM2 and NMF offer better results, despite some imperfections. The proposed methods behave similarly and are most likely to provide better results than the classical methods, in terms of estimation accuracy of the sources. Our methods maintain moderate numerical complexities if the number of observations is not too large. However, their complexities grow very fast as the number of observations increases. Fortunately, by incorporating a nonnegative compression preprocessing step, the overall complexity of the proposed SeNICA procedure is acceptable. These BSS applications demonstrate the validity and potential usefulness of the proposed algorithms. The experimental results have also shown that by fully exploiting the prior information of data, such as the nonnegativity of the mixing matrix and the statistical independence of the sources, the proposed methods achieve better estimation results.

Forthcoming work

The future work is twofold. The first part consists of improving the NJDC methods proposed in this PhD dissertation. The second part is to develop novel BSS methods which are more suitable for separating partially correlated sources.

Throughout the chapter 2, our effort is not only to develop five Jacobi-like NJDC algorithms but also to reduce their potential numerical complexities by combining them with classical JDC algorithms. Nevertheless, the proposed Jacobi-like NJDC algorithms are still more costly than the classical unconstrained JDC algorithms, even though the numerical complexities of our methods are lower than that of one of the existing NJDC algorithms, which is based on the LM optimization. The high complexity reported in chapters 3 and 4 seems the main limitation of our methods, and restricts their usage especially for high dimensional data with a large number of observations and of sources. Fortunately, besides these Jacobi-like methods, the proposed NJDC algorithm based on the ADMM optimization and the direct-fit JDC criterion requires a smaller computational load which is similar to that of classical algorithms, and maintains competitive performance in terms of estimation accuracy. It seems that the ADMM optimization provides a good compromise between performance and complexity especially for large scale data. Consequently, we would like to continue to derive novel NJDC methods based on the minimization of the indirect-fit JDC criteria by using the ADMM approach. Since the indirect-fit JDC criteria involve less variables than the direct-fit criterion, we can expect that the forthcoming algorithms could be more efficient than the existing one. It is noteworthy that we will face two difficulties. The first potential difficulty is to impose the nonnegativity on the joint transformation matrix A, meanwhile performing the ADMM optimization on its inverse A -1 . The second one is to avoid the trivial solution A -1 = 0 in the ADMM procedure. In the proposed Jacobi-like methods based on the indirect-fit criteria, the nonsingularity of A is guaranteed by the LU and QR decompositions, since the resulting matrix has a unit-determinant. In the ADMM method with the indirect-fit criteria, an additional nonsingularity constraint on A seems necessary. Moreover, further theoretical studies of the convergence of these approaches will be considered. On the other hand, recent studies on nonorthogonal Jacobi-like JDC algorithms showed that the convergence rate of the Jacobi-like procedure can be accelerated by considering the coupled LU and QR matrix factorizations [Maurandi et al., 2013, Maurandi and Moreau, 2014a[START_REF] Maurandi | A new fast Jacobi-like algorithm for non-orthogonal joint diagonalization of real-valued matrices based on a QR parameterization[END_REF], Maurandi et al., 2014]. We would like to develop new Jacobi-like NJDC algorithms by using such coupled factorizations in order to accelerate convergence and then to reduce the numerical complexity.

In chapter 4, we illustrate the usefulness of our SeNICA methods in the BSS framework. The success of SeNICA, as well as ICA and NICA, relies on the statistical independence, or at least on the quasi-statistical independence of the sources. However, in many practical problems, this assumption does not always hold. Hence, ICA, NICA and our SeNICA methods look for latent sources that are as statistically independent as possible. The degree of independence is measured through the use of mathematical tools from HOS, information theory and many others. Estimating the sources by minimizing the statistical dependence may not give satisfactory results, especially when the true source signals are known to be correlated [START_REF] Naanaa | Blind source separation of positive and partially correlated data[END_REF]. This situation appears in a number of real-world applications, such as hyperspectral unmixing [START_REF] Chan | A convex analysis-based minimum-volume enclosing simplex algorithm for hyperspectral unmixing[END_REF], Ma et al., 2014] and medical imaging [Wang et al., 2010]. Regarding the MRS data that we treated in this thesis, sometimes the source spectra can be highly correlated if the molecules share common structural features and therefore are difficult to separate. Consequently, an accurate separation of correlated or dependent sources still remains a challenging task. In [Wang et al., 2010], Wang et al. proposed a nonnegative Least-correlated Component Analysis (nLCA) method to separate nonnegative correlated sources by minimizing a joint correlation objective function among the estimated sources. Such an objective function assumes that the nonnegative sources are linearly independent rather than statistically independent, and exploits the convex geometry of the nonnegative mixtures of nonnegative sources. nLCA was shown to have a superior performance on real biomedical data over several existing benchmark methods, for instance, ICA and NMF. However, this method is constrained to work on a square mixing matrix and is sensitive to additional noise. We are working on a novel method which generalizes the nLCA algorithm to overdetermined mixtures of sources. Preliminary experimental results show that this novel method is more robust to additional noise and is also computationally efficient. Further studies on separating correlated MRS signals will be considered in order to systematically evaluate the performance of the new method.

where Q 1 = 2 t i,j /(1 + 2β j t i,j )βe T j , where β is defined in equation (2.85), and where β j is the j-th element of vector β. In the following, we compute the components of the terms 1 ○, 2 ○, 3 ○ and 4 ○ in equation (A.5). In 1 ○, T (i,j) (-t 2 i,j ) and T (i,j) (-t 2 i,j ) T transform the i-th row vector and i-th column vector of matrix C (k) as follows: In 2 ○, matrix Q 1 C (k) can be expressed as follows:

1 ○ =               i-th column | C (k) 1,1 C (k
Q 1 C (k) = 2 t i,j 1 + 2β j t i,j βe T j C (k) = 2 t i,j 1 + 2β j t i,j C (k) j,1 β, C (k)
j,2 β, . . . , C

j,N β (A.7)

T (i,j) (-t 2 i,j ) and T (i,j) (-t 2 i,j ) T transform the i-th row vector and i-th column vector of matrix Q 1 C (k) as follows:

2 ○ = 2 t i,j 1 + 2β j t i,j ×                  i-th column | C (k) j,1 β 1 C (k) j,
C (k) j,1 β N C (k) j,2 β N . . . -C (k) j,j β N t 2 i,j + C (k) j,i β N . . . C (k) j,N β N                  (A.8) Term 3
○ is the transpose of 2 ○. Its expression is therefore omitted. In 4 ○, matrix Q 1 C (k) Q T 1 can be expressed as follows:

Q 1 C (k) Q T 1 =
2 t i,j 1 + 2β j t i,j βe T j C (k) e j β T 2 t i,j 1 + 2β j t i,j = 4 C (k) j,j t 2 i,j

(1 + 2β j t i,j ) 2 ββ T (A.9)

T (i,j) (-t 2 i,j ) and T (i,j) (-t 2 i,j ) T transform the i-th row vector and i-th column vector of Then the minimum of (B.8) with respect to ∆θ can be reached at:

∆θ = - ∂ 2 Ψ(θ (it) i,j ) ∂(θ (it) i,j ) 2 -1
∂Ψ(θ (it) i,j ) ∂θ have the same sign. In this situation it is not necessary to perform the optimization since it will not decrease the criterion [START_REF] Ouedraogo | Non-negative independent component analysis algorithm based on 2D Givens rotations and a Newton optimization[END_REF]. The numerical complexity of estimating one θ i,j costs (9N + 1)J (it) θ i,j , where J (it) θ i,j denotes the number of consumed iterations. The process of all the rotation angles θ i,j is called a rotation sweep. The rotation sweep can be stopped when the relative decrease of the total sum of squares of negative components of V between two sweeps falls below a fixed small threshold.

B.2.2 Step 2: estimation of the elementary upper triangular matrices

After the rotation step, the elementary upper triangular matrix U (i,j) (u i,j ) transforms each pair of row vectors (i, j) of V as follows: W j,: W i,:

= 1 u i,j 0 1

V j,: V i,: (B.13)
where W i,: and W j,: are the i-th and j-th row vectors of matrix W respectively, and where u i,j is also called a shearing factor. From equation (B.13), the vector W i,: remains unchanged. The negativity measure criterion is then defined as follows:

Ψ(u i,j ) = 1 2 N n=1 W 2 j,n 1 W j,n <0 (B.14)
Similarly, the second order Taylor expansion Ψ T (u i,j ) of Ψ(u i,j ) around u

  source separation problem . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Independent component analysis . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Identifiability and preprocessing . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Non-Gaussianity, maximum likelihood and mutual information . . 3 1.2.3 ICA by tensor decomposition of cumulant arrays . . . . . . . . . . 6 1.3 From ICA to joint diagonalization by congruence . . . . . . . . . . . . . . 1.3.1 JDC problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Uniqueness and identifiability . . . . . . . . . . . . . . . . . . . . . 1.3.3 JDC algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Nonnegative matrix factorization . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Existence and uniqueness . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Basic NMF algorithms . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.3 Semi-NMF algorithms . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Nonnegative independent component analysis . . . . . . . . . . . . . . . . 1.5.1 Preprocessing and identifiablity . . . . . . . . . . . . . . . . . . . . 1.5.2 NICA algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Semi-nonnegative ICA and nonnegative JDC . . . . . . . . . . . . . . . . 1.6.1 The semi-nonnegative ICA problem . . . . . . . . . . . . . . . . . . 1.6.2 Nonnegative JDC algorithms . . . . . . . . . . . . . . . . . . . . . 1.6.3 Contribution and outline of thesis . . . . . . . . . . . . . . . . . . . 1.7 List of publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

•

  CPD: Canonical Polyadic Decomposition • BSS: Blind Source Separation • ICA: Independent Component Analysis • INDSCAL: INdividual Differences in SCALing • JDC: Joint Diagonalization by Congruence • NJDC: Nonnegative Joint Diagonalization by Congruence • MI: Mutual Information • ML: Maximum Likelihood • NICA: Nonnegative Independent Component Analysis • NMF: Nonnegative Matrix Factorization • PCA: Principal Component Analysis • SeNICA: Semi-Nonnegative Independent Component Analysis Statistics • SO: Second Order • TO: Third Order • FO: Fourth Order • HO: Higher Order • HOS: Higher Order Statistics • iid: independent and identically distributed • pdf: probability density function Methods • AC: Alternating Columns • ACDC: Alternating Columns and Diagonal Center • ADMM: Alternating Direction Method of Multipliers • ALM: Augmented Lagrangian Multiplier • ALS: Alternating Least Squares • ANLS: Alternating Nonnegative Least Squares • CLU: Coupled LU matrix factorization • DC: Diagonal Center • DIEM: DIagonalization using Equivalent Matrices • EJD: Exact Joint Diagonalization • ELS: Enhanced Line Search • EVD: EigenValue Decomposition • FFDIAG: Fast Frobenius DIAGonalization • GEVD: Generalized EigenValue Decomposition • HALS: Hierarchical Alternating Least Squares • InforMax: Information Maximization • JAD: Joint Approximate Diagonalization • JADE: Joint Approximate Diagonalization of Eigen-matrices • LM: Levenberg Marquardt • MU: Multiplicative Updates • NLS: Nonnegativity constrained Least-Square • NN-COMP: NonNegative COMPression • QDIAG: Quadratic DIAGonalization • SOBI: Second Order Blind Identification • STOTD: Simultaneous Third-Order Tensor Diagonalization • SVD: Singular Value Decomposition • UWEDGE: Uniformly Weighted Exhaustive Diagonalization with Gauss itErations Others • CPM: Completely Positive Matrix • ECG: ElectroCardioGraphy • EEG: ElectroEncephalography • flop: floating point operation • GB: GigaByte • GKLD: Generalized Kullback-Leibler Divergence • KKT: Karush-Kuhn-Tucker • MEG: MagnetoEncephaloGraphy • MRS: Magnetic Resonance Spectroscopy • NP: Non-deterministic Polynomial time • SED: Square of Euclidean Distance • SNR: Signal to Noise Ratio

•

  cum{.} denotes the cumulant function Usual matrices • 0 N : N -dimensional column vector of zeros • 0 N,P : (N × P )-dimensional matrix of zeros • I N : (N × N )-dimensional identity matrix Variables • N : number of observations • P : number of sources • M : number of sample points • K: number of matrices to be jointly diagonalized Norms • |.|: absolute value • . : Euclidean norm • . F : Frobenius norm IX List of Figures 3.1 The convergence in terms of minimizing the cost functions for estimating a (3 × 3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3) with three different SNR values. . 3.2 The convergence in terms of minimizing the cost functions for estimating a (5 × 3) non-square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (5 × 5) with three different SNR values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Influence of SNR for estimating a (3 × 3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3). 3.4 Influence of SNR for estimating a (5 × 3) non-square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (5 × 5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Influence of number K of the input symmetric matrices C (k) of size (3×3) for estimating a (3 × 3) square nonnegative joint transformation matrix A with SNR = 10 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Influence of number K of the input symmetric matrices C (k) of size (5×5) for estimating a (5 × 3) square nonnegative joint transformation matrix A with SNR = 10 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Influence of coherence of the column vectors of A for estimating a (3 × 3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3) with SNR = 10 dB. . . . . . . . . . . . . . . 3.8 Influence of coherence of the column vectors of A for estimating a (5 × 3) non-square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (5 × 5) with SNR = 10 dB. . . . . . . . . 3.9 Influence of modulus of uniqueness for estimating a (3 × 3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3) with SNR = 10 dB. . . . . . . . . . . . . . . . . . . . . 3.10 Influence of modulus of uniqueness for estimating a (5×3) non-square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (5 × 5) with SNR = 10 dB. . . . . . . . . . . . . . . . . . . . . X 4.1 An example of the BSS results of 2 simulated MRS metabolites by the proposed four methods without using the nonnegative compression step. The number of observations is set to N = 20 and the SNR value is fixed to 10 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Continuation of the first example of the BSS results of 2 simulated MRS metabolites by the proposed six methods with the nonnegative compression step. The number of observations is set to N = 20 and the SNR value is fixed to 10 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Continuation of the first example of the BSS results of 2 simulated MRS metabolites by five classical BSS methods. The number of observations is set to N = 20 and the SNR value is fixed to 10 dB. The red dot circles indicate the significant estimation errors in figures (b) to (f). . . . . . . . 4.4 Influence of the number of observations for estimating 2 MRS source metabolites with SNR = 10 dB. . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Influence of SNR for estimating 2 MRS source metabolites from 20 observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Three representations of the music data. The music waveform is contaminated by pink noises with a SNR value of 5 dB. . . . . . . . . . . . . . . . 4.7 Separated spectra obtained from six BSS methods from decomposing the music magnitude spectrogram. The red dot circles indicate the significant estimation errors in figures (d) to (f). . . . . . . . . . . . . . . . . . . . . . 4.8 Temporal weights obtained by six BSS methods by decomposing the music magnitude spectrogram. Each row in each subfigure represents the amplitude of the corresponding column vector of matrix A. . . . . . . . . 4.9 Restoration of the grayscale front and back sides of a real document by six BSS methods. The input images are obtained from [Tonazzini et al., 2007]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 An illustration of the proposed NN-COMP algorithm. (a) the input prewhitening matrix W ; (b) and (c) the proposed algorithm. (b) V is obtained by multiplying W by a Givens rotation matrix. (c) W is obtained by multiplying V by an elementary upper triangular matrix. The blue circles denote the nonnegative values and the red circles denote the negative values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . Moreover, they introduced a way to parameterize a nonnegative matrix A ∈ R N ×P + over the open set R N ×P by means of a square change of variable: A = B B, where B ∈ R N ×P has unconstrained components.

  52) namely a square change of variable and an exponential change of variable, respectively, where B ∈ R N ×P and where exp(•) denotes the element-wise exponential operator of the matrix argument. Then the nonnegativity constrained problem was transformed into a unconstrained problem. As previously reviewed in section 1.4.2, the square change of 1.6. Semi-nonnegative ICA and nonnegative JDC Chapter 1

  Now, let us return to the NJDC problem 1.8, where the nonnegativity constraint of the joint transformation matrix A ∈ R N ×N + is imposed through a square change variable A = B 2 . By means of the LU and QR matrix factorizations of B ∈ R N ×N , the jointly diagonalizable matrices C (k) = AD (k) A T can be rewritten in the following two equivalent forms, respectively:

  Figure 3.1: The convergence in terms of minimizing the cost functions for estimating a (3 × 3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3) with three different SNR values.

  Figure 3.2: The convergence in terms of minimizing the cost functions for estimating a (5 × 3) non-square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (5 × 5) with three different SNR values.

  Figure 3.3: Influence of SNR for estimating a (3 × 3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3).

  Figure 3.5: Influence of number K of the input symmetric matrices C (k) of size (3×3) for estimating a (3 × 3) square nonnegative joint transformation matrix A with SNR = 10 dB.

  Figure 3.6: Influence of number K of the input symmetric matrices C (k) of size (5×5) for estimating a (5 × 3) square nonnegative joint transformation matrix A with SNR = 10 dB.

  Figure 3.7: Influence of coherence of the column vectors of A for estimating a (3 × 3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3) with SNR = 10 dB.

  Figure 3.8: Influence of coherence of the column vectors of A for estimating a (5 × 3) non-square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (5 × 5) with SNR = 10 dB.

  Figure 3.9: Influence of modulus of uniqueness for estimating a (3×3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3) with SNR = 10 dB.

  Figure 3.10: Influence of modulus of uniqueness for estimating a (5 × 3) non-square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (5 × 5) with SNR = 10 dB.

3.

  Estimating the mixing matrix A ∈ R N ×P + or the compressed mixing matrix A + = WA ∈ R P ×P + by applying the proposed NJDC methods on the FO cumulant matrix slices of C (4) x or C (4)

  Estimation accuracy of the source vector s as a function of the number of observations Numerical complexity as a function of the number of observations

Figure 4 . 4 :

 44 Figure 4.4: Influence of the number of observations for estimating 2 MRS source metabolites with SNR = 10 dB.

  Figure 4.5: Influence of SNR for estimating 2 MRS source metabolites from 20 observations.

  Figure 4.7: Separated spectra obtained from six BSS methods from decomposing the music magnitude spectrogram. The red dot circles indicate the significant estimation errors in figures (d) to (f).

Figure 4 . 8 :

 48 Figure 4.8: Temporal weights obtained by six BSS methods by decomposing the music magnitude spectrogram. Each row in each subfigure represents the amplitude of the corresponding column vector of matrix A.

  4

  Figure 4.9: Restoration of the grayscale front and back sides of a real document by six BSS methods. The input images are obtained from[START_REF] Tonazzini | Fast correction of bleed-through distortion in grayscale documents by a blind source separation technique[END_REF].

  ) 2 = 0 means that all the elements of the row vectors V

  is the solution at the it-th iteration, ∆u = u i,j -u ) 2 are the first and second order derivatives of the cost function (B.14) with respect to u (it) i,j , respectively, which are given as follows:

Table 1 .

 1 2: The differences between the proposed six nonnegative JDC methods according to the cost functions, the optimization schemes and the assumption on A. When we use the square change variable A = B 2 , the LU and QR matrix factorizations are used to decompose B, leading to the first five Jacobi-like optimization methods. The sixth method based on ADMM optimization estimates A

	without change

  [START_REF] Wang | Nonnegative compression for semi-nonnegative independent component analysis[END_REF]. Nonnegative compression for semi-nonnegative independent component analysis. In SAM'14, Proceedings of the Eighth IEEE Sensor Array and Multichannel Signal Processing Workshop, pages 81-84, A Coruna,[Coloigner et al., 2014b] Coloigner, J., Fargeas, A.,Kachenoura, A., Wang, L., Senhadji, L., De Crevoisier, R., Acosta, O., and Albera, L. (2014b). A novel classification method for prediction of rectal bleeding in prostate cancer radiotherapy based on a semi-nonnegative ICA of 3D planned dose distributions. Accepted for publication in IEEE Journal of Biomedical and Health Informatics.[[START_REF] Wang | Parameterization of the inverse of the nonnegative matrix for nonnegative joint diagonalization[END_REF]] Wang, L., Albera, L.,[START_REF] Wang | Parameterization of the inverse of the nonnegative matrix for nonnegative joint diagonalization[END_REF]. Parameterization of the inverse of the nonnegative matrix for nonnegative joint diagonalization. In preparation for submission to IEEE Transactions on Signal Processing.

	Chapter 1
	Spain.
	3. International journals
	[Wang et al., 2013a] Wang, L., Albera, L., Kachenoura, A., Shu, H. Z., and
	Senhadji, L. (2013a). Nonnegative joint diagonalization by congruence based on
	LU matrix factorization. IEEE Signal Processing Letters, 20(8):807-810.

1. International conference without proceedings

[START_REF] Wang | Decomposition of semi-nonnegative semi-symmetric three-way tensors based on LU matrix factorization[END_REF]

] Wang, L.,

[START_REF] Wang | Decomposition of semi-nonnegative semi-symmetric three-way tensors based on LU matrix factorization[END_REF]

. Decomposition of semi-nonnegative semi-symmetric three-way tensors based on LU matrix factorization. In SIAM 2012 Conference on Applied Linear Algebra, Valencia, Spain.

2. International conferences with proceedings

[

Wang et al., 2013b

] Wang, L.,

Albera, L., Shu, H. Z., and Senhadji, L. (2013b)

. A new Jacobi-like nonnegative joint diagonalization by congruence. In EUSIPCO'13, Proceedings of the XXI European Signal Processing Conference, pages 1-5, Marrakech, Morocco.

[

[START_REF] Wang | CP decomposition of semi-nonnegative semi-symmetric tensors based on QR matrix factorization[END_REF]

] Wang, L., Albera, L.,

[START_REF] Wang | CP decomposition of semi-nonnegative semi-symmetric tensors based on QR matrix factorization[END_REF]

. CP decomposition of semi-nonnegative semi-symmetric tensors based on QR matrix factorization. In SAM'14, Proceedings of the Eighth IEEE Sensor Array and Multichannel Signal Processing Workshop, pages 449-452, A Coruna, Spain. (invited paper).

[

[START_REF] Wang | Nonnegative compression for semi-nonnegative independent component analysis[END_REF]

] Wang, L.,

Kachenoura, A., Albera, L., Karfoul, A., [Wang et al., 2014a

] Wang, L., Albera, L.,

Kachenoura, A., Shu, H. Z., and Senhadji, L. (2014a)

. Canonical polyadic decomposition of 3rd order semi-nonnegative semi-symmetric tensors using LU and QR matrix factorizations. Accepted for publication in EURASIP Journal on Advances in Signal Processing.

  By such a manipulation, most JDC algorithms based on criterion (1.27) can now estimate A directly. However, none of them can guarantee the nonnegativity of A. In order to impose the nonnegativity constraint on A, we resort to use a square change of variable A = B 2 . Then cost function (2.46) can be reformulated including the nonnegativity constraint as follows:

	Chapter 2	2.3. The iJDC + LU-1 and iJDC + QR-1 algorithms
		46)

  independent of the N -P vectors e p , therefore Å is nonsingular. Similarly, we can define the following three (N × N ) matrices:B = [B, e P +1 , . . . , e N ](2.120) L = [L, e P +1 , . . . , e N ]

	2.4. The JDC + LU-1 and JDC + LU-2 algorithms	Chapter 2

119)

Equation (2.119) suggests that we can transform the non-square matrix A into a square matrix Å by combining A with additional columns {e p } p ∈ {P + 1, P + 2, . . . , N }, without changing the structure of C

(k)

. In general, each column vector a p is linearly

  Table 3.1: Average numerical complexities of twelve methods in the convergence test for estimating a (3 × 3) square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (3 × 3) with various SNR values. Table 3.2: Average numerical complexities of six methods in the convergence test for estimating a (5 × 3) non-square nonnegative joint transformation matrix A from 5 input symmetric matrices C (k) of size (5 × 5) with various SNR values.

		0237 × 10 5	5.3056 × 10 4	5.0350 ×
	CLU	7.7080 × 10 3	1.5691 × 10 4	7.4604 ×
	FFDIAG	8.9549 × 10 3	1.2108 × 10 4 5.9718 × 10 3
	LUJ1D	3.2279 × 10 3	9.1461 × 10 3	2.2162 ×
	LUJ2D	7.8555 × 10 3	1.9839 × 10 4	4.5952 ×
	LM + sym	2.1544 × 10 6	5.4142 × 10 6	3.7247 ×
	iJDC + LU-1	1.8195 × 10 5	5.7355 × 10 4	1.2168 × 10 4
	iJDC + QR-1	1.7893 × 10 6	1.1665 × 10 6	1.2487 ×
	JDC + LU-1	4.5719 × 10 5	1.2280 × 10 6	7.6316 ×
	JDC + LU-2	1.9668 × 10 5	1.3533 × 10 5	6.8720 ×
	JDC + LU-3	1.6055 × 10 6	1.1272 × 10 6	7.9567 ×
	JDC + ADMM-3	1.0582 × 10 5	6.1208 × 10 4	6.9601 ×
		Computational complexity (flops)
		SNR = 30 dB SNR = 10 dB SNR = 0 dB
	ACDC	1.2473 × 10 5	9.2324 × 10 4	5.5712 × 10 4
	LM + sym	2.0315 × 10 6	9.2810 × 10 6	6.6496 ×
	JDC + LU-1	1.0388 × 10 6	4.0737 × 10 6	4.3547 ×
	JDC + LU-2	4.8495 × 10 5	6.8402 × 10 5	5.1805 ×
	JDC + LU-3	3.2191 × 10 6	4.1906 × 10 6	3.7078 ×
	JDC + ADMM-3	8.2344 × 10 4	7.8530 × 10 4	9.1460 × 10 4

  .3. Document restoration Table 4.3: Numerical complexities of six BSS methods in the experiment of document restoration.

		Chapter 4
		Numerical complexity (flops)
	CoM2	3.3017 × 10 7
	NICA	7.7118 × 10 8
	NMF	1.0613 × 10 11
	Semi-NMF	1.4150 × 10 8
	JDC + LU-2 -ICA	2.1279 × 10 7
	JDC + ADMM-3 -ICA	2.1272 × 10 7

Nous avons développé six nouvelles méthodes de DCC nonnégative, qui ont été soigneusement présentées dans ce chapitre. Les six algorithmes proposés peuvent être répartis dans deux groupes selon la stratégie d'optimisation. Le premier groupe d'algorithmes comprend cinq méthodes semi-algébriques, reposant sur une méthode de type Jacobi. Dans ce groupe, la contrainte de nonnégativité a été assurée par un changement de variable carré, permettant ainsi de nous ramener à un problème d'optimisation sans contrainte. L'idée générale de notre approche de type Jacobi est i) de factoriser la racine carrée de la matrice de passage recherchée sous forme d'un produit de matrices élémentaires, chacun n'étant défini que par un seul paramètre, puis ii) d'estimer ces matrices élémentaires l'une après l'autre dans un ordre spécifique. La première méthode utilise la factorisation LU afin de minimiser un critère des moindres carrés directs. Cependant, elle présente une grande complexité numérique. Afin de réduire la complexité numérique, nous avons reformulé le critère des moindres carrés indirects et nous avons considéré les factorisations LU et QR. Une telle reformulation a nécessitée l'utilisation d'une étape de prétraitement sur les matrices à diagonaliser, conduisant à deux algorithmes de calcul efficace. Néanmoins, la validité de ces deux algorithmes est fondée sur certaines hypothèses d'inverisibilité. Par conséquent, les quatrième et cinquième algorithmes proposés ont permis de relâcher ces hypothèses. Ces deux méthodes reposent sur une optimisation des deux critères des moindres carrés indirects basée sur la factorisation LU, sans aucune étape de prétaitement. Le deuxième groupe compte un seul algorithme, qui utilise la méthode des directions alternées. Cet algorithme exploite le fait que le problème de DCC nonnégative peut être reformulé comme un problème des moindres carrés directs sous contraintes de symétrie et de nonnégativité. Un tel algorithme a été dérivé en minimisant séquentiellement le Lagrangien augmenté par rapport aux variables et aux multiplicateurs. En

Chapter 2

The JDC +

LU-1 and JDC + LU-2 algorithms

For both algorithms, the initialization of C (k) in equations (2.86) requires 2N 3 K flops.

Regarding the JDC + LU-1 algorithm, in order to compute one parameter t i,j in one Jacobilike iteration, computing the coefficient matrix Q E in the cost function (2.98) involves (23N 2 -4N -3)K + N 3 + N 2 + N flops. Calculating the roots of the rational function (2.98) costs 729 flops. The computation load of A (new) , B (new) and C (k,new) , with k ∈ {1, • • • , K}, is 4N + (2N 2 + 5N + 2)K flops. Otherwise, if the parameter t i,j is computed without the nonnegativity constraint by means of equation (2.126), its calculation load is 2(N -1)K flops. Accordingly, the computation cost of A (new) and C (k,new) is 2N K + N flops. The row balancing scheme needs 3N 2 K + N P flops. Thus, the magnitude of global numerical complexity of JDC + LU-1 consisting of estimating all the P (N -1) parameters is between O(N 2 P K) and O(N 4 P + N 3 P K), depending on the number of unconstrained optimizations is performed. Regarding the JDC + LU-2 algorithm, in each Jacobi-like iteration, the estimation of t i,j by minimizing the cost function (2.115) requires (5N 3 + 32N 2 + 11N -5)K + 2N 3 + 5N 2 + 4N + 1331 flops. Computing A (new) , B (new) and C (k,new) involves 4N + (2N 2 + 5N + 2)K flops. Nevertheless, if t i,j is estimated by minimizing the unconstrained cost function (2.127), it costs 8K + 64 flops. In this case, updating A (new) and C (k,new) costs 2N K + N flops. Consequently, considering all the P (N -1) parameters, the magnitude of the total numerical complexities of JDC + LU-2 is between O(N 2 P K) and O(N 4 P K). The full expressions of the complexities of the JDC + LU-1 and JDC + LU-2 algorithms are given in table 2.1.

The JDC + ADMM-3 algorithm

Computing A (it+1) 1

according to equation (2.136) requires (N 2 P + N P + 2P 2 )K + N P 2 + P 3 + N P + P flops. The complexity of A (it+1) 2

is the same as that of A (it+1) 1

. Updating U (it+1) by means of equation (2.138) involves 3N P flops. According to equation (2.139), the calculation of K diagonal matrices D (k,it+1) needs N 2 P K + N 2 P 2 + N 2 P + N P 2 + P 3 +P 2 flops. Finally, the computation of the multipliers Π (it+1) 1 and Π (it+1) 2 both require N P flops, following equations (2.140) and (2.141), respectively. Therefore the magnitude of the total numerical complexities of one iteration of JDC + ADMM-3 is O(N 2 P K), while its full expression is given in table 2.1.

From table 2.1, we can observe that in order to estimate a square joint transformation matrix A ∈ R N ×N + , the classical JDC methods ACDC, CLU, FFDIAG, LUJ1D and LUJ2D achieve lower numerical complexities per sweep/iteration than most of the nonnegativity constrained methods. The magnitude of the complexities of these classical JDC methods is O(N 3 K). Among the proposed methods, the JDC + ADMM-3 algorithm gives the lowest complexity, which is also of magnitude O(N 3 K). The complexities of the JDC + LU-3 , iJDC + LU-1 and iJDC + QR-1 algorithms are at most one magnitude higher than that of JDC + ADMM-3 , which are of magnitude O(N 4 K) in the worst case when all the parameters are computed by the nonnegativity constrained optimizations. The JDC + LU-1 and JDC + LU-2 algorithms require a large amount of computation, whose magnitudes of complexities are O(N 5 + N 4 K) and O(N 5 K) in the worst case, respectively. Nevertheless, their complexities are still much lower that of the Levenberg Marquardt method, namely LM + sym , when K is much larger than N . The complexity of LM + sym is of magni-Chapter 2

2.6. Numerical complexity analysis tude O(N 3 K 3 + N 6 ) due to the heavy computation of the modification of the Hessian. A similar observation can be made from table 2.1 when a non-square joint transformation matrix A ∈ R N ×P + with N > P is encountered in a JDC problem.

Chapter 3

3.4. Influence of number of input matrices algorithm is illustrated in figure 3.4(c). One can see that generally LM + sym has the heaviest computation load. JDC + LU-1 and JDC + LU-3 are more efficient than LM + sym under most SNR levels. ACDC requires the lowest complexity when SNR is less than 8 dB, while JDC + ADMM-3 remains the most computationally efficient when SNR exceeds 8 dB. The numerical complexity of JDC + LU-2 is between those of ACDC and JDC + LU-1 . Thus, JDC + LU-2 and JDC + ADMM-3 provide the best accuracy/speed compromise.

Influence of number of input matrices

In this section, we study the influence of the number K of the input matrices C (k) on the performance of the considered algorithms. In ICA, the value of K corresponds for instance to the number of matrix slices derived from a FO cumulant array.

Test on square matrices

In this experiment, the dimensions of the square joint transformation matrix A ∈ R N ×P + are set to N = 3 and P = 3. The SNR value is fixed to 10 dB. We repeat the experiment with K ranging from 3 to 63 with a step of 2. Figure 3.5(a) shows the average curves of the estimation error α(A, A) of all the algorithms as a function of K. We can observe that the performance of all the methods progresses as K increases, and then practically stabilizes for high values of K. It indicates that after some point (e.g. K ≥ 30), the additional information brought by an increase of K does not further improve the results.

Regarding the four classical Jacobi-like methods, CLU, FFDIAG, LUJ1D and LUJ2D, their performance is quite similar. ACDC outperforms the aforementioned four algorithms. All the nonnegativity constrained methods give better results than the five classical unconstrained algorithms. The proposed JDC + LU-1 and JDC + LU-2 algorithms maintain competitive advantages over the other methods through all the K values. iJDC + QR-1 and JDC + ADMM-3 are also among the most competitive methods, and they achieve better results than LM + sym under most K values. iJDC + LU-1 shows similar performance to LM + sym . The performance of JDC + LU-3 is between that of ACDC and LM + sym , and it is also comparable. The average numerical complexities of all the algorithms in this experiment are presented in figure 3.5(b). It shows that the numerical complexities of all the methods progress quasi-linearly as K increases except that of LM + sym , which grows quadratically in K. Now let's give a detailed description. The complexities of the classical Jacobi-like methods, CLU, FFDIAG, LUJ1D and LUJ2D, remain the lowest, especially FFDIAG. The proposed iJDC + LU-1 , JDC + LU-2 and JDC + ADMM-3 algorithms cost a bit more flops than ACDC. The complexities of iJDC + QR-1 , JDC + LU-1 and JDC + LU-3 are almost one order of magnitude higher than that of ACDC. Nevertheless, they are still much more efficient than LM + sym , especially for a large K value. Therefore, the JDC + LU-2 and JDC + ADMM-3 methods seem to be the most effective algorithms compared to the other methods in terms of both estimation accuracy and numerical complexity.

Test on non-square matrices

In this experiment, the dimensions of the non-square joint transformation matrix A ∈ R N ×P + are set to N = 5 and P = 3, and the SNR value is fixed to 10 dB. Again, we repeat the experiment with K ranging from 3 to 63 with a step of 2. The average curves Chapter 3

Influence of coherence of matrix A

where each column vector a p ∈ R N + is randomly drawn from a uniform distribution over [0, 1], and where β is in the range of [0, 1]. A small value of β indicates that the matrix A is close to be singular.

Test on square matrices

In this simulation, the dimensions of the square joint transformation matrix A ∈ R N ×P + are set to N = 3 and P = 3, the number of the input matrices C (k) is fixed to K = 5, and the SNR value is set to 10 dB. We repeat the experiment with β ranging from 0 to 1 with a step of 0.025. The average curves of the estimation error α(A, A) of the twelve algorithms are depicted in figure 3.7(a). One can notice that in general all the JDC algorithms are less effective in terms of the estimation accuracy when two column vectors of A are close to collinear. The performance of most JDC methods increases quasi-linearly as β increases. The four classical Jacobi-like methods CLU, FFDIAG, LUJ1D and LUJ2D yield less accurate estimation results. ACDC exhibits an improved performance over the aforementioned methods. The nonnegativity constrained algorithms give much better results than the classical unconstrained methods. Three methods based on the direct-fit JDC criterion (1.29), namely LM + sym , JDC + LU-3 and JDC + ADMM-3 , show similar performance. The three proposed Jacobi-like NJDC algorithms, namely iJDC + LU-1 , iJDC +

QR-1

and JDC + LU-2 , outperform the previous three methods. The performance of JDC + LU-1 is quite stable and remains the best through all the considered β values. Figure 3.7(b) depicts the average numerical complexities of all the algorithms in this experiment. We can observe that the complexity evolution is quite stable as β changes. LM + sym is the most computationally expensive, followed by JDC + LU-3 , JDC + LU-1 and iJDC + QR-1 , consecutively. iJDC + LU-1 and JDC + ADMM-3 cost a bit more or less flops than ACDC, while JDC + LU-2 is somewhat more expensive than ACDC. The classical Jacobi-like methods maintain the lowest numerical complexities. Hence, JDC + LU-2 gives the best accuracy/complexity compromise in this experiment, and JDC + LU-1 is also more efficient than LM + sym .

Test on non-square matrices

In this experiment, the dimensions of the non-square joint transformation matrix A ∈ R N ×P + are set to N = 5 and P = 3, and the number of the input matrices C (k) is fixed to K = 5. The SNR value is set to 10 dB. We repeat the experiment with β ranging from 0 to 1 with a step of 0.025. The average curves of the estimation error α(A, A) of the six considered algorithms as a function of SNR are presented in figure 3.8(a). The obtained results show that the nonnegativity constrained methods outperform the classical ACDC algorithm. The direct-fit JDC criterion (1.29) based methods, LM + sym , JDC + LU-3 and JDC + ADMM-3 give similar results, where the proposed JDC + ADMM-3 algorithm is slightly better than the other two algorithms. JDC + LU-1 gives quite stable performance as β increase, and remains advantage when β is less than 0.4. JDC + LU-2 generates the most accurate results when β exceeds 0.5. The average numerical complexities evolutions of all the six methods as a function of β are illustrated in figure 3.8(b). We can see that ACDC and JDC + ADMM-3 are the most computationally efficient, followed by JDC + LU-2 , JDC + LU-1 and JDC + LU-3 , consecutively. LM + sym requires the heaviest workload. Therefore, the proposed JDC + LU-2 and JDC + ADMM-3 algorithms offer the best performance/complexity compromise in this experiment.

Test on non-square matrices

In this experiment, the dimensions of the non-square joint transformation matrix A ∈ R N ×P + are set to N = 5 and P = 3, and the number of the input matrices C (k) is fixed to K = 5. The SNR value is fixed to 10 dB. We repeat the experiment with the angle ψ ranging from 0 to π/2 with a step of π/60. Figure 3.10(a) displays the average curves of the estimation error α(A, A) of the six selected algorithms as a function of ψ. One can observe that the nonnegativity constrained algorithms surpass the classical ACDC algorithm, except JDC + LU-1 when ψ > 5π/12. Three methods based on the direct-fit JDC criterion (1.29), called LM + sym , JDC + LU-3 and JDC + ADMM-3 behave similarly, where the proposed JDC + ADMM-3 algorithm offers better performance and gives the best estimation result when ψ > 2π/5. The proposed two Jacobi-like algorithms JDC + LU-1 and JDC + LU-2 , particularly JDC + LU-1 , achieve large performance gains when ψ is less than 3π/10. The evolutions of the average numerical complexities of the six methods as a function of ψ are depicted in figure 3.10(b). It shows that ACDC and JDC + ADMM-3 keep the lowest computational complexities, while ACDC is the most efficient when ψ < 2π/10. The complexity of JDC + LU-2 is almost one order of magnitude higher than that of JDC + ADMM-3 . JDC + LU-1 and JDC + LU-3 need the heaviest workload among the proposed methods. Nevertheless, they are still more efficient than LM + sym , since the complexities of the former are half of the latter. Therefore, the proposed JDC + LU-2 algorithm gives the best performance/complexity compromise in this experiment.

Concluding remark

In this chapter, the performance of the proposed six NJDC algorithms was evaluated with simulated set of jointly diagonalizable matrices. Five classical nonorthogonal JDC methods without using the nonnegativity constraint, including ACDC [Yeredor, 2002], CLU [Maurandi et al., 2013], FFDIAG [Ziehe et al., 2004], LUJ1D [Afsari, 2006], LUJ2D [Afsari, 2006], and one NJDC method LM + sym [Coloigner et al., 2014c], were tested as reference methods. The performance was assessed in terms of the matrix estimation accuracy and the numerical complexity. Then, various scenarios were considered, aiming at testing: i) the convergence property, ii) the influence of additive noise, which can stand for model errors, iii) the impact of the number of input matrices, iv) the effect of coherence of the column vectors of the joint transformation matrix, and v) the influence of the modulus of uniqueness. In each scenario, the performance was extensively studied through Monte Carlo experiments for both estimating a square joint transformation matrix and a non-square one. The obtained results showed that the proposed algorithms offer better estimation accuracy in difficult contexts by means of exploiting the nonnegativity a priori. Comparing to an existing nonnegative method based on the Levenberg Marquardt optimization [Coloigner et al., 2014c], the proposed methods reduce the total numerical complexity by up to three orders of magnitude. For a difficult problem with either a lower value of SNR, a small number of available input matrices, a high level of coherence presented in the joint transformation matrix or a large value of modulus of uniqueness, generally the proposed JDC + LU-2 algorithm provides the best performance/complexity compromise. JDC + LU-1 also gives the best estimation results in many difficult contexts, where a highly coherent joint transformation matrix is encountered or a large Chapter 4

Reconstructing the M realizations of the source vector by means of s

, where denotes the pseudo-inverse operator. The complexity of this step is 2N P 2 + N P M + P 3 or P 2 M + P 3 by using the above two equations, respectively.

The first step is indispensable for the iJDC + LU-1 -ICA and iJDC + QR-1 -ICA algorithms when a non-square mixing matrix is encountered, since they only work on a square matrix. Such a compression step is optional for the JDC + LU-1 -ICA, JDC + LU-2 -ICA, JDC + LU-3 -ICA and JDC + ADMM-3 -ICA methods. The numerical complexities of the proposed six SeNICA methods include the complexities of the four steps. Their full expressions are straightforward, thus omitted here for brevity. The behaviors of the proposed SeNICA methods are compared with five state-of-the-art BSS algorithms, namely two efficient ICA methods CoM2 [Comon, 1994] and SOBI [Belouchrani et al., 1997], the NICA method with a line search along the geodesic [Plumbley, 2003], the NMF method based on alternating nonnegativity least squares [Kim and Park, 2008], and the semi-NMF method based on one-sided multiplicative update [Ding et al., 2010]. A comprehensive analysis of the numerical complexities of the classical BSS algorithms is provided in [START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF], Albera et al., 2010, Albera et al., 2012]. Their results are not repeated here for the sake of conciseness.

Separation of simulated MRS signal

MRS is a powerful non-invasive analytical technique for analyzing the chemical content of Magnetic Resonance (MR)-visible nuclei and therefore enjoys particular advantages for assessing metabolism. The chemical property of each nucleus determines the frequency at which it appears in the MR spectrum, giving rise to peaks corresponding to specific metabolites [START_REF] Befroy | Magnetic resonance spectroscopy studies of human metabolism[END_REF]. Therefore, the MRS observation spectra can be modeled as the mixture of the spectrum of each constitutional source metabolite. More specifically, it can be written as the noisy linear instantaneous mixing model described in equation (1.50) as in the SeNICA model 1.6, where x ∈ R N is the MRS observation vector, s ∈ R P is the source vector representing the statistically quasi-independent source metabolites, ε ∈ R N is the instrumental noise vector, and A ∈ R N ×P + is the nonnegative mixing matrix containing the concentrations of the source metabolites. In this context, SNR is defined as follows:

Many non-parametric BSS methods were developed in order to solve the MRS separation problem, such as the ICA approach [START_REF] Nakai | Application of independent component analysis to magnetic resonance imaging for enhancing the contrast of gray and white matter[END_REF], Szabo de Edelenyi et al., 2005] and the NMF approach [Sajda et al., 2004, Sun and[START_REF] Sun | [END_REF].

In this experiment, two simulated MRS source metabolites s 1 and s 2 , namely the Choline and Myo-inositol (see figure 4.1(b)), are generated by the sum of Lorentzian and Gaussian functions, where the location and scale parameters are fixed for the purposed of realistic representation [Moussaoui, 2005]. Each of the sources contains 10 3 samples. The observation vector x is generated according to equation (1.50) (see figure 4.1(a) for illustration of two observations with a SNR of 10 dB). The components of the (N × 2) mixing matrix A are randomly drawn from a uniform distribution. The additive noise ε is modeled as a zero-mean unit-variance Gaussian vector. The performance is assessed by means of the estimation error α(s T , s T ) between the true source s and its estimate s, and the numerical complexity, where the distance measure α is defined in equation (3.2). Furthermore, all the results reported in the section are averaged over 200 independent Monte Carlo trials. 1994], NICA [Plumbley, 2003], NMF [Kim and Park, 2008], semi-NMF [Ding et al., 2010] and SOBI [Belouchrani et al., 1997] Table 4.1: The theoretical fundamental frequencies of four musical notes, namely C 5 , A b 4 , F 4 and D b 4 , on an ideal piano [Jorgensen, 1991].

4.6(b) and 4.6(c), respectively. Two of the proposed methods are chosen to solve this problem, namely JDC + LU-2 -ICA and JDC + ADMM-3 -ICA, since they provide a good compromise between estimation accuracy and numerical complexity. It is noteworthy that computing the FO cumulant array of the magnitude spectrogram X ∈ R 251×1379 + requires at least (251 4 × 8)/2 30 ≈ 29.6 GB memories, which is beyond the capacity of ordinary computers. Therefore, for the proposed methods, we propose to use NN-COMP as a preprocessing step. The proposed JDC + LU-2 -ICA and JDC + ADMM-3 -ICA incorporating NN-COMP are compared with CoM2 [Comon, 1994], NICA [Plumbley, 2003], NMF [Kim and Park, 2008] and semi-NMF [Ding et al., 2010]. Figure 4.7 displays the separated spectra obtained by the considered six methods. It can be seen that the proposed NN-COMP + JDC , respectively. The fifth component seems to capture the transient events corresponding to the attacks and releases of the notes, more precisely, the sound produced by the hammer hitting the string and by the release of the sustain pedal. The sixth one appears to consist of the pink noise spectrum. 4 and C 5 , and note F 4 appears twice in the separated spectra. Figure 4.8 illustrates the temporal weight matrix A obtained by the six methods, where each row in each subfigure represents the amplitude of the corresponding column vector of A. The locations and widths of the peaks appearing in each column vector of A characterize the the onset times and durations of one note. Therefore, the correctness of the result can be verified by checking if the temporal weights coincide with the original music sheet. For the proposed NN-COMP + JDC + LU-2 -ICA and NN-COMP + JDC + ADMM-3 -ICA methods, as well as CoM2, the locations of the peaks in the first four rows of figures 4.8(b) to 4.8(d) well match the original music sheet in figure 4.8(a), while the widths of these peaks seem identical since all the notes have the same duration. It is worth noting that in figure 4.8(d) the temporal weights of note C 5 obtained by CoM2 contain a large amount of negative values, which could bring some difficulties for the physical interpretation. The nonnegativity constraint of the weight matrix A helps the results of the proposed methods giving a more realistic physical meaning. Concerning NICA and NMF, the estimated weights for note C 5 are obviously incorrect, as illustrated in figures 4.8(e) and 4.8(f). Semi-NMF gives the worst result, since all the peaks in figure 4.8(g) are covered by noises. Additionally, the numerical complexities of all the six methods are listed in table 4.2. CoM2 has the lowest complexity. The complexities of the proposed NN-COMP + JDC + LU-2 -ICA and NN-COMP + JDC + ADMM-3 -ICA methods are moderately higher than that of CoM2, and between those of semi-NMF and NMF. Therefore, in this experimental context, the proposed SeNICA algorithms considering both the statistical independence of the source spectra and the nonnegativity of the temporal weights improve the separation result and also provide a good accuracy/complexity tradeoff.

Proofs of propositions

A.1 Proof of proposition 2.1

Multiplying B by an elementary triangular matrix T (i,j) (t i,j ) on the right side performs the elementary transformation of the j-th column vector of B:

where b i , i ∈ {1, 2, . . . , N } is the i-th column vector of B. By routine computation, the Hadamard square of equation (A.1) is:

A.2 Proof of proposition 2.2

It is straightforward to show that the elements of the j-th column vector except the (j, j)-th entry of the term 1 ○ in equation (2.51) can be expressed by

n,j with n ∈ {1, . . . , N } and n = j, and those elements of the term 2 ○ in (2.51) are equal to c (k,1) n u i,j with n ∈ {1, . . . , N } and n = j. The sum of these elements directly leads to equation (2.52). The terms 3 ○ and 4 ○ in equation (2.51) do not need to be considered, since they do not affect the off-diagonal elements in the j-th column vector.

A.3 Proof of proposition 2.3

Multiplying B by a Givens rotation matrix Q (i,j) (θ i,j ) on the right side performs the elementary transformation of the i-th and j-th column vectors of B:

3) where b i , i ∈ {1, 2, . . . , N } is the i-th column vector of B. By routine computation, the Hadamard square of equation (A.3) is:

i -e T j ).

A.4 Proof of proposition 2.4

It is straightforward to show that the (i, j)-th entry of the term 1 ○ in (2.61) can be expressed by sin 2 (θ i,j ) cos 2 (θ i,j )( C

), and that of the term 4 ○ in (2.61) is -sin 2 (2θ i,j ) c (k,3) . Then proposition 2.4 can be proved.

A.5 Proofs of propositions 2.6 and 2.7

Let's expand equation (2.86) as follows:

10) The sum of equations (A.6), (A.8), (A.10) and the transpose of (A.8) will directly lead to propositions 2.6 and 2.7.

A.6 Proof of proposition 2.8

According to proposition 2.1, the differences between A (new) and A are their j-th column vectors. By the definition of Khatri-Rao product [Cichocki et al., 2009, Chapter 1], the following equality exists:

where a j is the j-th column vector of A. Hence, the differences between A (new) A (new) and A A are their j-th column vectors. The j-th column vector of A (new) A (new) can be expressed as follows:

where {g i } 4 i=1 , G and ζ i,j are defined in proposition 2.8. The (N 2 × 1) column vector vec ( a j a T j ) is the j-th column vector of A A. Then proposition 2.8 holds.

A.7 Proof of proposition 2.9

We can rearrange the ADMM update formulas in equations (2.136) to (2.141) into the following form:

(A.13) A.7. Proof of proposition 2.9

)d (k,it) (A. 16) .18) where c (k) = vec(C (k) ), and where d (k,it) , d (k,it+1) ∈ R P are column vectors containing the diagonal elements of the diagonal matrices D (k,it) and D (k,it+1) , respectively. The assumption Z (it+1) -Z (it) → 0 implies that the left and right hand sides of equations (A.13) to (A.18) go to zero. Let it go to infinity. We have

1 ), ... , and so on for each variable, where the second terms

1 , ... , vanishes asymptotically. Consequently, we obtain:

where the terms α(U (it+1) -A (it+1) 1

) and β(U (it+1) -A (it+1) 2

) have been eliminated from equations (A.19) and (A.20), respectively, by invoking equation (A.23) and (A.24). Obviously, the first five equations (2.149) to (2.153) in the KKT conditions of problem (2.129) are satisfied at any limit point:

From equations (A.23) and (A.24), we have

Therefore, the following condition can be proved:

Hence, proposition 2.9 holds.

Nonnegative compression method for SeNICA

B.1 Problem formulation

In some practical SeNICA problems defined in equation (1.50), the dimension of the observation space x must be reduced. The classical dimension compression procedure, such as the spatial prewhitening, compresses the N -dimensional vector x into a vector x of dimension P N . The estimate of the rank P is determined by the number of eigenvalues of the covariance matrix of x not exceedingly close to zero. The compressed observation vector x is expressed as follows [Comon, 1994, Belouchrani et al., 1997]:

where W ∈ R P ×N is called a prewhitening matrix. In order to compute W , we obtain the square root of the covariance matrix, denoted by Γ ∈ R N ×P , such that ΓΓ T = E{xx T } assuming x being centered. Then W = Γ where denotes the pseudo-inverse operator [Belouchrani et al., 1997]. However, the prewhitening matrix W breaks the nonnegativity property of the compressed mixing matrix A = W A. Thus the proposed NJDC algorithms cannot be applied to estimate A. In this section, we introduce a new Non-Negative COMPression method, namely NN-COMP, which guarantees the nonnegativity of the compressed mixing matrix A. The NN-COMP algorithm computes a linear transformation matrix Φ ∈ R P ×P such that the resulting compression matrix W = ΦW has nonnegative components. Consequently, now the compressed observation vector x can expressed as follows:

where the compressed mixing matrix A + = WA ∈ R P ×P is a nonnegative square matrix.

Then the proposed NJDC algorithms can be used to compute the compressed matrix A + from the HO cumulant array of x. Once A + is estimated, the original matrix A is obtained as follows:

3)

It should be noted that the source s can be recovered by s = A -1 + x. Therefore generally it is not necessary to compute the original mixing matrix A. Such a novel compression method can reduce the numerical complexity of the proposed NJDC algorithms. It is particularly important for the iJDC + LU-1 and iJDC + QR-1 algorithms, since they are constrained to work on a square matrix A. Now, let's define the nonnegative compression problem as follows:

Problem B.1. Given a prewhitening matrix W ∈ R P ×N of an N -dimensional random vector process x, find a sequence of Givens rotation matrices and elementary upper triangular matrices, such that their product:

have nonnegative components. Then the transformation matrix Φ ∈ R P ×P is defined as the product of all these Givens rotation matrices and elementary upper triangular matrices in the proper order. A Givens rotation matrix searches for a rotation angle that makes the outputs matrix V as nonnegative as possible. However, sometimes it still remains some negative values near the quadrant boundaries (figure B.1(b)). That is because the row vectors of W are neither well-grounded nor statistically independent [Plumbley, 2003]. We propose to use an elementary upper triangular matrix which projects the remaining negative values of V into the nonnegative quadrant (figure B.1(c)). We thus obtain a nonnegative compression matrix W and the (P × P ) compressed mixing matrix A + = W A, which preserves the nonnegativity property. Now the challenge is how to compute the Givens rotation and the elementary upper triangular matrices.

B.2. Algorithm derivation

B.2 Algorithm derivation

A new two-step NN-COMP method (figure B.1) is presented in this section.

B.2.1 Step 1: estimation of the Givens rotation matrices

The Givens rotation matrix R (i,j) (θ i,j ) transforms each pair of row vectors (i, j) of W as follows:

V i,:

V j,: = cos(θ i,j ) sin(θ i,j ) -sin(θ i,j ) cos(θ i,j ) W i,: W j,:

(B.5)

where V i,: (V j,: ) and W i,: (W j,: ) are the i-th (j-th) row vectors of matrices V and W , respectively. Now, let us consider the following negativity measure criterion defined in [Plumbley, 2003, Ouedraogo et al., 2010]:

where V i,n and V j,n are the (i, n)-th and (j, n)-th components of V , respectively, and where 1 α<0 is a Heaviside-step-like function defined as follows:

The purpose is to find an angle, θ i,j , which minimizes the total sum of squares of negative components of V . The global optimum of Ψ(θ i,j ) is difficult to obtain analytically due to the existence of the Heaviside-step-like function 1 α<0 . We propose to compute θ i,j

iteratively by a Newton's method. For a given iteration (it), let us consider the second order Taylor expansion, Ψ T (θ i,j ), of Ψ(θ i,j ) around θ (it) i,j :

where θ

(it) i,j is the solution at the it-th iteration, ∆θ = θ i,j -θ

i,j ) 2 are the first and second order derivatives of the cost function (B.6) with respect to θ (it) i,j , respectively, which are given as follows [START_REF] Ouedraogo | Non-negative independent component analysis algorithm based on 2D Givens rotations and a Newton optimization[END_REF]:

and where

Then the minimum of (B.15) with respect to ∆u can be reached at:

i,j ) 2 = 0, it means that all the elements of the row vector W (it) j,: have the same sign and that equation (B.19) does not need to be computed. The numerical complexity of estimating one u i,j costs (4N + 1)J (it) u i,j , where J (it) u i,j denotes the number of used iterations. The process of all the shearing factors u i,j is called a shearing sweep. The shearing sweep can be stopped when the relative decrease of the total sum of squares of negative components of W between two sweeps falls below a predefined positive threshold.

In practice, the new NN-COMP algorithm contains some rotation sweeps followed by several shearing sweeps, in order to ensure the convergence. The overall numerical complexity of NN-COMP is given as follows:

where N , P and M are the numbers of the observations, of the sources and of the sample points, respectively, where J (it) θ i,j and J (it) u i,j are the numbers of inner iterations for estimating one θ i,j and one u i,j , respectively, and where J Compute Θ (i,j) (ξi,j)) by equation (2.12) Initialization: k,0) and it = 0 3:

The AC phase 6:

for j = 1 to N do Estimate L 7:

for i = j + 1 to N do 8:

Compute i,j by minimizing (2.44)

9:

ς = sign( aj + i,j ai) for i = 1 to N do Estimate U 20:

Compute ui,j by minimizing (2.44)

22:

ς = sign( aj + ui,j ai) end for 32:

The DC phase 33:

Compute D (k,new) by (2.41)

35:

end for 36:

end while

38:

Return:

) -1 and it = 0

for j = 1 to N do Estimate L 6:

for i = j + 1 to N do 7:

Compute i,j by (2.72) new) 31:

end for 32:

end for 33:

if (it mod q) == 0 then Row balancing of C (k) 34:

Perform row balancing scheme

) -1 and it = 0

Compute θi,j by minimizing (2.73)

8:

ςi = sign(cos (θi,j) ai -sin (θi,j) aj), ςj = sign(sin (θi,j) ai + cos (θi,j) aj)

9:

if |(cos (θi,j) ai -sin (θi,j) aj)| == ςi(cos (θi,j) ai -sin (θi,j) aj) and 10:

|(sin (θi,j) ai + cos (θi,j) aj)| == ςj(sin (θi,j) ai + cos (θi,j) aj) then 11: new) 19:

end for 20:

end for 21:

Compute ui,j by (2.72)

24:

ς = sign( aj + ui,j ai) new) 33:

end for 34:

end for 35:

if (it mod q) == 0 then Row balancing of C (k) 36:

Perform row balancing scheme for i = j + 1 to N do 7:

Compute i,j by (2.126)

Compute ui,j by (2.126)

22:

ς = sign( aj + ui,j ai) if (it mod q) == 0 then Row balancing of C (k) 34:

Perform row balancing scheme for i = j + 1 to N do 7:

Compute i,j by minimizing (2.127)

8:

ς = sign( aj + i,j ai) if (it mod q) == 0 then Row balancing of C (k) 34:

Perform row balancing scheme Initialization:

2 , U (0) , {D (k,0) }, Π Compute U (it+1) by equation (2.138)

Compute each D (k,it+1) by equation (2.139)

Compute Π Initialization: compute the prewhitening matrix Γ such that ΓΓ T = E{xx T } by EVD;