Canonical polyadic decomposition of 3rd order semi-nonnegative semi-symmetric tensors using LU and QR matrix factorizations

Semi-symmetric three-way arrays are essential tools in Blind Source Separation (BSS) particularly in Independent Component Analysis (ICA). These arrays can be built by resorting to higher order statistics of the data. The Canonical Polyadic (CP) decomposition of such semi-symmetric three-way arrays allows us to identify the so-called mixing matrix, which contains the information about the intensities of some latent source signals present in the observation channels. In addition, in many applications, such as the Magnetic Resonance Spectroscopy (MRS), the columns of the mixing matrix are viewed as relative concentrations of the spectra of the chemical components. Therefore, the two loading matrices of the three-way array, which are equal to the mixing matrix, are nonnegative. Most existing CP algorithms handle the symmetry and the nonnegativity separately. Up to now, very few of them consider both the semi-nonnegativity and the semi-symmetry structure of the three-way array. Nevertheless, like all the methods based on line search, trust region strategies and alternating optimization, they appear to be dependent on initialization, requiring in practice a multi-initialization procedure. In order to overcome this drawback, we propose two new methods, called JD + LU and JD + QR , to solve the problem of CP decomposition of semi-nonnegative semi-symmetric three-way arrays. Firstly, we rewrite the constrained optimization problem as an unconstrained one. In fact, the nonnegativity constraint of the two symmetric modes is ensured by means of a square change of variable. Secondly, a Jacobi-like optimization procedure is adopted because of its good convergence property. More precisely, the two new methods use LU and QR matrix factorizations, respectively, which consist in formulating high-dimensional optimization problems into several sequential polynomial and rational subproblems. By using both LU and QR matrix factorizations we aim at studying the inﬂuence of the used matrix factorization. Numerical experiments on simulated arrays emphasize the advantages of the proposed methods especially the one based on LU factorization, in the presence of high-variance model error and of degeneracies such as bottlenecks. A BSS application on MRS data conﬁrms the validity and improvement of the proposed methods.


Introduction
Higher Order (HO) arrays, commonly called tensors, play an important role in numerous applications, such as chemometrics [1], telecommunications [2], and biomed-ical signal processing [3].They can be seen as HO extensions of vectors (1-way arrays) and matrices (2-way arrays).In many practical situations the available data measurements cannot be arranged into a tensor form directly, that is to say, the observation diversity is insufficient either in time or frequency.However, if the latent data satisfies the statistical independence assumption, which is reasonable in many applications, meaningful HO arrays can be built by resorting to HO Statistics (HOS) of the data [4].In this instance, the HO arrays are partially symmetric or Hermitian due to the special algebraic structure of the basic HOS, such as moments and cumulants.In Independent Component Analysis (ICA), the latent physical phenomena which are assumed to be statistically independent can be revealed by decomposing the HO array into factors.There exists several ways to decompose a given HO array, such as the Tucker model [5,6].Among the existing reliable HO array decomposition models, the Canonical Polyadic (CP) decomposition model have attracted much attention.Indeed its uniqueness properties can be ensured under the sufficient conditions established by Kruskal [7].In addition unlike the HO Singular Value Decomposition (HOSVD) [6], CP model does not impose any orthogonality constraint on its factors.
Theoretically, a polyadic decomposition exactly fits an array by a sum of rankone terms [8].A CP decomposition is defined as a polyadic decomposition with a minimal number of rank-one terms which are needed to exactly fit a given HO array.Currently the CP decomposition is gaining importance in several applications, for example, in exploratory data analysis [9], sensor array processing [10], telecommunications [11,12], ICA [13] and in multiple-input multiple-output radar systems [14].A multitude of methods were developed to compute the CP decomposition.They include the iterative Alternating Least Squares (ALS) procedure [15], which gains popularity due to its simplicity of implementation and low numerical complexity.Uschmajew proved the local convergence property of ALS under some conditions [16].However, this convergence can be slow.Therefore an Enhanced Line Search (ELS) procedure was proposed by Rajih et al. [17] to cope with the slow convergence problem of ALS.Other approaches were also proposed, such as the conjugate gradient algorithm [18] and joint eigenvalue decomposition based algorithms [19,20], to cite a few.Some HO arrays enjoy certain properties, such as i ) symmetry and ii ) nonnegativity, which can not be simply handled by the aforementioned general CP decomposition methods.Therefore, special CP models become more and more important.
The first special form of the CP model for 3-way arrays that are symmetric in two modes, brings the concept of INdividuals Differences in SCALing (INDSCAL) analysis [21].On one hand, INDSCAL analysis has been studied as a way of multiple factor analysis [22] with applications to chemometrics, psychology and marketing research.On the other hand, in the domain of signal processing, and more particularly in Blind Source Separation (BSS), the INDSCAL analysis is widely known as the Joint Diagonalization of a set of matrices by Congruence (JDC).During the past two decades, many successful JDC methods have been proposed, such as Yeredor's Alternating Columns and Diagonal Center (ACDC) algorithm [23], the Joint Approximate Diagonalization (JAD) algorithm proposed by Cardoso and Souloumiac [24], the Fast Frobenius DIAGonalization (FFDIAG) algorithm proposed by Ziehe et al. [25], Afsari's LUJ1D algorithm [26], and many others [27][28][29][30][31][32][33].
A recent survey of JDC can be found in [34].The second special form of CP model is defined when all the factors in the CP decomposition are constrained to be nonnegative, commonly known as Nonnegative Tensor Factorization (NTF).NTF can be regarded as the extension of Nonnegative Matrix Factorization (NMF) [35] to higher orders.In many applications, the physical properties are inherently nonnegative, such as chemistry [1] and fluorescence spectroscopy [36,37].In those applications the results are only meaningful if the nonnegativity constraint is satisfied.Various methods for computing NTF and also NMF can be found in [38,39].
So far, the CP model with both the symmetry and nonnegativity constraints has not received much attention.Coloigner et al. proposed a family of algorithms based on line search and trust region strategies [40].Wang et al. developed an alternating minimization scheme [41].Those methods appear to depend on initialization, and therefore in practice require a multi-initialization procedure, leading to an increase of numerical complexity.In this paper, we propose to fit the CP model of a 3-way array by imposing both the semi-nonnegativity and the semi-symmetry constraints.More precisely, we impose a nonnegativity constraint on the two symmetric modes of the INDSCAL model, which leads to the semi-nonnegative INDSCAL model or equivalently the CP decomposition of semi-nonnegative semi-symmetric 3-way arrays.Such a model is often encountered in ICA problems where a nonnegative mixing matrix is frequently considered.For example, in Magnetic Resonance Spectroscopy (MRS), the columns of the mixing matrix represent the positive concentrations of the source metabolites.Then the 3-way array built by stacking the matrix slices of a cumulant array is both nonnegative and symmetric in two modes.In such a case, the semi-nonnegative INDSCAL problem is equivalent to the JDC problem subject to a nonnegativity constraint on the joint transformation matrix.We propose two new algorithms to solve the semi-nonnegative INDSCAL problem, called JD + LU and JD + QR .Firstly, we rewrite the constrained optimization problem as an unconstrained one.Actually the nonnegativity constraint is ensured by means of a square change of variable.Secondly, we propose two Jacobi-like approaches using LU and QR matrix factorizations, respectively, which consist in formulating highdimensional optimization problems into several sequential polynomial and rational subproblems.By using both LU and QR matrix factorizations we aim at studying the influence of the used matrix factorization.Numerical experiments highlight the advantages of the proposed methods especially JD + LU , in the case of dealing with high-variance model error and with degeneracies such as bottlenecks.A BSS application on MRS signals confirms the validity and improvement of the proposed methods.A part of this work has been recently presented at the 8th IEEE Sensor Array and Multichannel signal processing workshop [42].
The rest of the paper is organized as follows.After the presentation of some notations, the second section introduces some basic definitions of the multilinear algebra then gives the semi-nonnegative INDSCAL problem formulation.In section three, we describe the proposed algorithms in details and we also provide an analysis of the numerical complexities.Section four shows the computer simulation results.Finally we conclude the paper.

Notations
The following notations are used throughout this paper.R N1×N2× The (i, j)-th entry of a matrix A is symbolized by A i,j .Sometimes the MATLAB column/row notation is adopted to indicate submatrices of a given matrix or subarrays of a HO array.Also a i denotes the i-th column vector of matrix A. denotes the Hadamard product (element-wise product), and A 2 = A A. denotes the Khatri-Rao product.A denotes the pseudo inverse of A. The superscripts, −1 , T and −T stand for the inverse, the transpose and the inverse after transpose operators, respectively.The (N × N ) identity matrix is denoted by I N .0 N stands for N -dimensional vectors of zeros.|a| denotes the absolute value of a.A F and det(A) stand for the Frobenius norm and determinant of matrix A, respectively.diag(A) returns a matrix comprising only the diagonal elements of A. Diag(b) is the diagonal matrix whose diagonal elements are given by the vector b.off(A) vanishes the diagonal components of the input matrix A. vec(A) reshapes a matrix A into a column vector by stacking its columns vertically.

Definitions and problem formulation
Now we introduce some basic definitions in multilinear algebra which are necessary for the problem formulation.

Definition 1
The outer product C = u (1) • u (2) • u (3) of three vectors u (i) ∈ R Ni (1 ≤ i ≤ 3) is a 3-way array of R N1×N2×N3 whose elements are defined by Definition 2 Each 3-way array C expressed as the outer product of three vectors is a rank-1 3-way array.
More generally, the rank of a 3-way array is defined as follows: Definition 3 The rank of an array C ∈ R N1×N2×N3 , denoted by rk(C), is the minimal number of rank-1 arrays belonging to R N1×N2×N3 that yield C in a linear combination.
Despite the similarity between the definition of the tensor rank and its matrix counterpart, the rank of a three-way array may exceed its dimensions [4].
Definition 4 A 3-way array slice is a 2-dimensional section (fragment) of a 3-way array, obtained by fixing one of the three indices [38].
For example, the k-th frontal slice of a 3-way array C can be denoted by C :,:,k using MATLAB notation and sometimes it is also denoted by C (k) .
The low-rank INDSCAL model of a 3-way array is defined as follows: Definition 5 For a given P , corresponding to the number of rank-1 terms, the INDSCAL model of a 3-way array C ∈ R N ×N ×K can be expressed as: where the 3-way array V represents the model residual.
If and only if the residual V is a null tensor, we have an exact INDSCAL decomposition.
An exact INDSCAL decomposition is considered to be essentially unique when it is only subject to scale and permutation indeterminacies.It means that an INDSCAL decomposition is insensitive to a scaling of the three vectors a p , a p and d p provided that the product of the 3 scale numbers is equal to 1, and an arbitrary permutation of the rank-1 terms.A necessary and sufficient uniqueness condition for the INDSCAL model was established by Afsari [43].
The INDSCAL model can also be described by using the frontal slices of C: where D (k) is a diagonal matrix whose diagonal contains the elements of the k-th row of D, and V (k) = V :,:,k .
In this paper we propose to fit the INDSCAL model of 3-way arrays, while imposing nonnegativity constraints on both equal loading matrices A. It will be referred to as the semi-nonnegative INDSCAL model, as follows: The semi-nonnegative INDSCAL problem is equivalent to the JDC problem subject to the nonnegativity constraint on the joint transformation matrix.In this paper, we mainly focus on the case of square nonnegative joint transformation matrix, for which N = P .The case of N > P will be discussed briefly in the next section.Therefore the problem that we tackle in this paper is defined as follows: by minimizing the residual term V (k) in a least-squares sense, subject to A having nonnegative components.

JDC cost functions
If the residual array V is a realization of a Gaussian random array, it is logical to fit the INDSCAL model by the following Direct Least Square (DLS) criterion [23,44]: and to minimize (4) with respect to A and D. Note that, in the field of ICA, only the loading matrix A is of interest since it corresponds to the mixing matrix of several latent source signals.The minimization of (4) with respect to D, when A is fixed, was given by Yeredor in [23]: When A is orthogonal, we can replace Then the extra parameter D can be eliminated and the minimization of ( 4) is equivalent to minimizing the following InDirect Least Square (IDLS) criterion [45,46]: In some cases such as in ICA, the orthogonality assumption of A can be satisfied by using a spatial whitening procedure [47].However, it is known that the whitening procedure may introduce additional errors.Therefore many algorithms propose to relax the orthogonality constraint by introducing the following cost function [25,31]: Frequently the minimization of criterion ( 7) is performed on a matrix Z def = A −1 instead of A for simplicity and Z is called the joint diagonalizer.To use this criterion, the matrix A (or Z) should be properly constrained in order to avoid the trivial zero solution and/or degenerate solutions [34].
Besides the criterions (4) and ( 7), Afsari [26] presented a new cost function, which is invariant to column scaling of A. Pham proposed an information theoretic criterion [48], which requires each matrix C (k) to be positive definite.Tichavský and Yeredor gave a special weighted least square criterion [49].

Problem reformulation
Existing semi-nonnegative INDSCAL algorithms are based on the minimization of the cost function (4) [40,41].They are able to achieve a better estimation of A than ACDC when the data satisfies the semi-nonnegative INDSCAL model at the cost of a higher computational complexity.We propose to use criterion (7) based on elementary factorizations of A due to the fast convergence property of this kind of procedures.Generally, it is quite difficult to impose the nonnegativity constraint on A, while computing its inverse A −1 by minimizing (7).Let us consider the structure of C = [[A, A, D]] with the following assumptions: 1 Then each frontal slice of C is nonsingular and its inverse can be expressed as follows: We use C (k,−1) to denote (C (k) ) −1 for simplicity.Equation (8) shows that C (k,−1) also preserves the jointly diagonalizable structure.Furthermore, instead of A −1 , A serves as the joint diagonalizer.Then A can be estimated by minimizing the following modified criterion based on (7): By such a manipulation, most algorithms based on criterion (7) 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 which was introduced by Chu et al. [50] for NMF, next adopted by Royer et al. for NTF [37] and by Coloigner et al. for semi-nonnegative INDSCAL [40]: where B ∈ R N ×N .Then problem 2 can be reformulated as follows: , find the square nonnegative loading matrix A = B 2 such that B minimizes the following cost function: LU and QR parameterizations of B In order to minimize (11), one may consider a gradient-like approach.However, the performance of this kind of method is sensitive to the initial guess and to the search step size.In addition, the calculation of gradient of (11) with respect to B is computationally expensive due to the existence of Hadamard product.Other algorithms, using Jacobi-like procedures [25,26,31], parameterize A as a product of several special elementary matrices and estimate each elementary matrix successively.We propose to follow such a minimization scheme.Now let us recall the following definitions and lemmas: Definition 6 A unit upper (or lower) triangular matrix is an upper (or lower, respectively) triangular matrix whose main diagonal elements are equal to 1.
Definition 7 An elementary upper (or lower) triangular matrix with parameters {i, j, u i,j } and i < j is a unit upper (or lower, respectively) triangular matrix whose non-diagonal elements are zeros except the (i, j)-th entry, which is equal to u i,j .
U (i,j) (u i,j ) with 1 ≤ i < j ≤ N denotes an elementary upper triangular matrix: Similarly, L (i,j) ( i,j ) with 1 ≤ j < i ≤ N corresponds to an elementary lower triangular matrix.
Definition 8 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.
) with 1 ≤ i < j ≤ N indicates the corresponding Givens rotation matrix: Lemma 1 Any (N × N ) unit lower triangular matrix L whose (i, j)-th component is i,j (i > j) can be factorized as the following product of N (N − 1)/2 elementary lower triangular matrices [51, Chapter 3]: where the two sets of indices J 1 and I 1 (j) are defined by J 1 = {1, 2, . . ., N } and I 1 (j) = {j + 1, j + 2, . . ., N } for the sake of convenience.Similarly, any (N × N ) unit upper triangular matrix U whose (i, j)-th component is equal to u i,j (i < j) can be factorized as a product of elementary upper triangular matrices as follows: where I 2 and J 2 (i) are two sets of indices, defined by Lemma 2 Any (N × N ) orthonormal matrix Q can be factorized as the following product of N (N − 1)/2 Givens rotation matrices [52,Chapter 14]: where I 2 and J 2 (i) are defined in lemma 1.
For any nonsingular matrix B ∈ R N ×N , the LU matrix factorization decomposes it as B = LU ΛΠ, where L ∈ R N ×N is a unit lower triangular matrix, U ∈ R N ×N is a unit upper triangular matrix, Λ ∈ R N ×N is a diagonal matrix and Π ∈ R N ×N is a permutation matrix.B also admits the QR matrix factorization as B = QR Λ, where Q ∈ R N ×N is an orthonormal matrix, R ∈ R N ×N is a unit upper triangular matrix, and Λ ∈ R N ×N is a diagonal matrix.Due to the indeterminacies of the JDC problem, the global minimum of ( 11), say B, can be expressed as B = LU and B = QR without loss of generality.Moreover, by incorporating lemma 1 and lemma 2, we obtain the two following elementary factorizations of B: As a consequence, the minimization of (11) with respect to B is converted to the estimate of N (N − 1) parameters: i,j and u i,j for the LU decomposition (17), or θ i,j and u i,j for the QR decomposition (18).Instead of simultaneously computing the N (N − 1) parameters, we propose two Jacobi-like procedures which perform N (N − 1) sequential optimizations.This yields two new algorithms: i ) the first algorithm based on (17), named JD + LU , estimates each i,j and u i,j successively, and ii ) the second one based on (18), called JD + QR , estimates each θ i,j and u i,j sequentially.Now, the difficulty is how to estimate four kinds of parameter, namely L (i,j) ( i,j ) and U (i,j) (u i,j ) for JD + LU , and Q (i,j) (θ i,j ) and U (i,j) (u i,j ) for JD + QR .Two points should be noted here: i ) 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 (i < j for U (i,j) (u i,j ) and j < i for L (i,j) ( i,j )); ii ) for both JD + LU and JD + QR algorithms, the procedure of estimating U (i,j) (u i,j ) is identical.Consequently, the principal problem is reduced to estimating two kinds of parameters, namely U (i,j) (u i,j ) and Q (i,j) (θ i,j ).
Minimization with respect to the elementary upper triangular matrix U (i,j) (u i,j ) In this section, we minimize (11) with respect to U (i,j) (u i,j ) with 1 ≤ i < j ≤ N .Let Ã and B denote the current estimate of A and B before estimating the parameter u i,j , respectively.Let Ã(new) and B(new) stand for Ã and B updated by U (i,j) (u i,j ), respectively.Furthermore, the update of B is defined as follows: In order to compute the parameter u i,j , a typical way is to minimize the criterion (11) with respect to u i,j by replacing matrix B by B(new) .For the sake of convenience, we denote J(u i,j ) instead of J( B(new) ).Then J(u i,j ) can be expressed as follows: The expression of the Hadamard square of the update B(new) is shown in the following proposition: )) 2 can be expressed as a function of u i,j as follows: where bi and bj denote the i-th and j-th columns of B, respectively, and e j is the j-th column of the identity matrix I N .
Inserting (21) into the cost function (20), we have: where constant column vector, a (1 × N ) constant row vector and a constant scalar, respectively.The term 1 ○ in ( 22) transforms the j-th column and the j-th 22) is a zero matrix except its j-th column containing non-zero elements, while the term 3 ○ contains non-zero entries only on its j-th row.And the term 4 ○ is a zero matrix except its (j, j)-th component being non-zero.In addition, Hence (22) shows that only the j-th column and j-th row of C(k,new) involve the parameter u i,j , while the other elements remain constant.Therefore, the minimization of the cost function ( 20) is equivalent to minimizing the sum of the squares of the j- The required elements of C(k,new) can be expressed by the following proposition.

Proposition 2
The elements of the j-th column except the (j, j)-th entry of C(k,new) is a second degree polynomial function in u i,j as follows, for every value n different of j: 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 straightforward.Indeed, we can show that the elements of the j-th column except the (j, j)-th entry of the term 1 ○ in ( 22) can be expressed by n,j with 1 ≤ n ≤ N and n = j, and those elements of the term 2 ○ in ( 22) are equal to c(k,1) n u i,j with 1 ≤ n ≤ N and n = j.The sum of these elements directly leads to (23).The terms 3 ○ and 4 ○ do not need to be considered, since they do not affect the off-diagonal elements in the j-th column.Proposition 2 shows that the minimization of the cost function ( 20) can be expressed in the following compact matrix form: where is a ((N − 1) × 3) matrix defined as follows: the first column contains the i-th 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.u i,j = [u 2 i,j , u i,j , 1] T is a 3-dimensional parameter vector.Equation (24) shows that J(u i,j ) is a fourth degree polynomial function.The global minimum u i,j can be obtained by computing the roots of its derivative and selecting the one yielding the smallest value of (24).Once the optimal u i,j is computed, B(new) is updated by (19) and the joint diagonalizer Ã(new) is updated by computing ( B(new) ) 2 .Then the same procedure is repeated to compute the next u i,j with another (i, j) index.
The minimization of (11) with respect to the elementary lower triangular matrix L (i,j) ( i,j ) with 1 ≤ j < i ≤ N can be computed in the same way.Proposition 2 is also valid for the parameter i,j when 1 ≤ j < i ≤ N .The detailed derivation is omitted here.The processing of all the N (N −1) parameters u i,j and i,j is called a LU sweep.In addition, for estimating L (i,j) ( i,j ), the (i, j) index obeys the following order: (2, 1), (3, 1), . . ., (N, 1),(3, 2), (4, 2), . . ., (N, 2), . . ., Regarding U (i,j) (u i,j ), the (i, j) index varies according to the following sequence: The proposed JD + LU algorithm is comprised of several LU sweeps.
Minimization with respect to the Given rotation matrix Q (i,j) (θ i,j ) Now we minimize (11) with respect to Q (i,j) (θ i,j ) with 1 ≤ i < j ≤ N .By abuse of notation, in this section we continue to use Ã and B to denote the current estimate of A and B, respectively, before estimating the parameter θ i,j .Also let Ã(new) and B(new) stand for Ã and B updated by Q (i,j) (θ i,j ), respectively.The update of B is defined as follows: Similarly, for computing the parameter θ i,j , we can minimize the criterion (11) with respect to θ i,j by replacing matrix B by B(new) .We denote J(θ i,j ) instead of J( B(new) ) for convenience purpose.Then J(θ i,j ) can be expressed as follows: The Hadamard square of the update B(new) now can be rewritten as shown in the following proposition.
)) 2 can be written as a function of θ i,j as follows: where bi and bj denote the i-th and j-th columns of B, respectively, and e i and e j are the i-th and j-th columns of the identity matrix I N , respectively.
Inserting (29) into the cost function (28), we obtain: where and c(k,3) = ( bi bj ) T C (k,−1) ( bi bj ) are a (N × N ) constant matrix, a (N × 1) constant column vector, a (1 × N ) constant row vector and a constant scalar, respectively.The term 1 ○ in (30) transforms the i-th and j-th columns and the i-th and j-th rows of C(k) .The term 2 ○ in ( 30) is a zero matrix except its i-th and j-th columns containing non-zero elements, while the term 3 ○ contains non-zero entries only on its i-th and j-th rows.And the term 4 ○ is a zero matrix except its (i, i)-th, (j, j)-th, (i, j)-th and (j, i)-th components being non-zero.
○ is a (N × N ) symmetric matrix.Hence (30) shows that only the i-th and j-th columns and the i-th and j-th rows 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 global cost function (28), such a simplified minimization scheme is commonly adopted in many algorithms, such as [20,31].We denote this local minimization by J(θ i,j ).The (i, j)- Proposition 4 The (i, j)-th entry of C(k,new) can be expressed as a function of θ i,j as follows: where 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.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.
It is straightforward to show that the (i, j)-th entry of the term 1 ○ in (30) can be expressed by sin 2 ), and that of the term 4 ○ is − sin 2 (2θ i,j )c (k,3) .Then proposition 4 can be proved.
In order to simplify the notation of (31), we resort to the Weierstrass change of variable: t i,j = tan(θ i,j ).Then we obtain: By substituting (32) into (31), we obtain an alternative expression of the (i, j)-th entry of C(k,new) which is described in the following proposition.Then the minimization of J(θ i,j ) transforms to J(t i,j ).
Proposition 5 the (i, j)-th entry of C(k,new) can be expressed by a rational function of t i,j as follows: where f Equation ( 33) easily shows that the sum of the squares of the (i, j)-th entries of the K matrices C(k,new) , is a rational function in t i,j , namely J(t i,j ), where the degrees of the numerator and the denominator are 8 and 8, respectively.J(t i,j ) can be expressed in the following compact matrix form: where T is a 5-dimensional vector, and τ i,j is a 5-dimensional parameter vector defined as follows, 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 J(t i,j ).Once t i,j is obtained, θ i,j can be computed from the inverse tangent function θ i,j = arctan(t i,j ).It is noteworthy that the found θ i,j cannot guarantee to decrease the actual cost function (28).If θ i,j leads to an increase of (28), we reset θ i,j = 0. Otherwise, B(new) is updated as described in (27) and the joint diagonalizer Ã(new) is updated by computing ( B(new) ) 2 .The same procedure will be repeated to compute θ i,j with the next (i, j) index.The order of the (i, j) indices is defined in equation (26).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 JD + QR algorithm.

Both of the JD +
LU and JD + QR algorithms can be stopped when the value of cost function (11) 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.

Practical issues
In practice, we observe that if each frontal slice of the 3-way array C is almost exactly jointly diagonalizable due to a high Signal to Noise Ratio (SNR), the classical non-constrained JDC methods can also give a nonnegative A with high probability.In this situation the explicit nonnegativity constraint could be unnecessary and could increase the computational burden.Therefore, we propose to relax the nonnegativity constraint by directly decomposing A into elementary LU and QR forms, respectively, instead of using the decompositions of B as follows: where the index sets I 1 (j), J 1 , I 2 and J 2 (i) are defined in lemma 1.By inserting ( 36) and ( 37) into the cost function ( 9), the ways of estimating the two sets of parameters { i,j , u i,j }, and {θ i,j , u i,j } are identical to those of Afsari's LUJ1D and QRJ1D methods [26], respectively.Therefore, in practice, in order to give an automatically SNR-adaptive method, for JD + LU , in each Jacobi-like iteration, we suggest to compute u i,j by LUJ1D first.If all the elements in the j-th column of ÃU (i,j) (u i,j ) have the same sign ε, the update Ã(new) = ε ÃU (i,j) (u i,j ) is adopted.
Otherwise, u i,j is computed by minimizing (20) and Ã(new) is updated by computing (21).Each i,j is computed similarly.Furthermore, the proposed JD + QR and QRJ1D are combined in the same manner.
Afsari reported in [26] that if the rows of matrices balanced in their norms, the computation of the parameter 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 Λ and Ã(new) by Ã(new) = ÃΛ using a diagonal matrix Λ ∈ R N ×N + , whose diagonal elements are defined as follows: where C(k) n,: denotes the n-th row of C(k) .
In ICA, when a non-square matrix A ∈ R N ×P + with N > P is encountered, the invertibility assumption of the frontal slices C (k) does not hold.In this situation, we can compress A by means of a nonnegative matrix W + ∈ R P ×N + such that the resulting matrix Ā = W + A is a nonnegative square matrix.Then the JD + LU and JD + QR algorithms can be used to compute the compressed loading matrix Ā. W + can be computed by using the NonNegative COMPression algorithm (NN-COMP) that we proposed in [53].More precisely, given a realization of an observation vector, we obtain the square root of the covariance matrix, denoted by Υ ∈ R N ×P .The classical prewhitening matrix is computed by W = Υ ∈ R P ×N where denotes the pseudo inverse operator [47].Then the NN-COMP algorithm computes a linear transformation matrix Ψ ∈ R P ×P such that W + = ΨW has nonnegative components.Once Ā is estimated, the original matrix A is obtained as follows: It should be noted that generally A does not need to be computed in such an ICA problem, since the sources can be estimated directly by means of Ā.

Numerical complexity
The numerical complexities of JD + LU and JD + QR are 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 loading matrix A ∈ R N ×N + from a 3-way array C ∈ R N ×N ×K , is considered; which does not affect the order of magnitude of the numerical complexity.
For both algorithms, the inverses C (k,−1) (k ∈ {1, • • • , K}) of the frontal slices of Ãini requires 2N 3 K flops, and at each sweep, the calculation of parameters u i,j needs N (N − 1)(5N 2 + 12N − 8)K/2 flops.In addition, in the case of the JD + LU algorithm, the calculation cost of Ã(new) , B(new) and flops, and the numerical complexity of computing the parameters i,j is equal to that of u i,j .Regarding the JD + QR algorithm, for each sweep, the complexity of calculating the parameters θ i,j is equal to N (N − 1)(5N  [41] is also based on a square change of variable and LU matrix factorization.It minimizes the cost function (4) with respect to A and D alternately, leading to a higher numerical complexity.By means of the reformulation of the cost function, the proposed methods avoid the estimation of D, therefore achieve a lower complexity compared to ACDC + LU .The explicit expressions of the overall complexity of JD + LU , JD + QR and ACDC + LU [41], as well as those of four classical JDC algorithms, namely ACDC [23], FFDIAG [25], LUJ1D [26] and QRJ1D [26], are listed in table 1.One can notice that numerical complexities of the proposed JD + LU and JD + QR methods are at most one order of magnitude higher than those of the four JDC algorithms, and still lower than that of ACDC + LU .Moreover, JD + LU is less computationally expensive than JD + QR .

Simulation results
This section is twofold.In the first part, the performance of the proposed JD + LU and JD + QR algorithms is evaluated with simulated semi-nonnegative semi-symmetric 3way arrays C. Several experiments are designed to study the convergence property, the influence of SNR, the impact of the third dimension K of C, the effect of the coherence of the loading matrix D and the influence of the condition number of the diagonal matrices D (k) .We also evaluate the proposed methods for estimating a non-square matrix A. The proposed algorithms are compared with four classical nonorthogonal JDC methods, namely ACDC [23], FFDIAG [25], LUJ1D [26], QRJ1D [26], and the nonnegative JDC method ACDC + LU [41].In the second part, the source separation ability of the proposed algorithms is studied through a BSS application.In this context, the JD + LU and JD + QR are used to jointly diagonalize several matrix slices of the fourth order cumulant array [40] of the observations and compared with several classical ICA [47,54,55] and NMF [56] methods.
where σ N is a scalar controlling the noise level.Then the SNR is defined by SNR = −20 log 10 (σ N ).ii) Initialization: in each Monte Carlo trial, all the algorithms are initialized by a same random matrix whose components obey the uniform distribution over [0, 1].iii) Afsari's row balancing scheme: the LUJ1D, QRJ1D, JD + LU and JD + QR algorithms perform the row balancing scheme once per each run of five sweeps.iv) Stopping criterion: all the algorithms stop either when the relative error of the corresponding criterion between two successive sweeps is less than 10 −5 or when the number of sweeps exceeds 200.A sweep of ACDC includes a full AC phase and a DC phase.v) Performance measurement: the performance is measured by means of the error between the true loading matrix A and the estimate Ã, the numerical complexity, and the CPU time.We define the following scale-invariant and permutation-invariant distance [40]: where a n and ãn are the n-th column of A and the n -th column of Ã, respectively.I 2 n is defined recursively by In addition, d(a n , ãn ) is defined as the pseudo-distance between two vectors [13]: The criterion ( 41) is an upper bound of the optimal permutation-invariant criterion.It avoids the burdensome computation of all the permutations.A small value of ( 41) means a good performance in the sense that Ã is close to vi) Test environment: 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 GB memory.Moreover, we repeat all the experiments with 500 Monte Carlo trials.

Convergence
In this experiment, the convergences of the JD + LU and JD + QR algorithms are compared to those of ACDC, FFDIAG, LUJ1D, QRJ1D and ACDC + LU QR , it reduces the cost function (11) to the values relatively higher than those achieved by JD + LU , and converges in about 50 sweeps whatever the SNR is.It seems that FFDIAG, LUJ1D and QRJ1D achieve the fastest convergence rate.It should be noted that while an algorithm may converge to a point in which the value of the cost function is close to zero, such a point could be a local minimum far from the desire matrix A as shown in figure 2. The top picture in figure 2 shows the convergence curves measured in terms of the estimating error α(A, Ã) as a function of sweeps when SNR = 25 dB.It shows that the solutions of FFDIAG, LUJ1D and QRJ1D are still far from optimum.ACDC and ACDC + LU give better estimations of A than the previous three methods.The best results are achieved by the proposed JD + LU and JD + QR methods.The middle picture in figure 2 displays the convergence curves when SNR = 10 dB.It can be observed that ACDC converges to a local minimum which is not the global one, and that the performance of the proposed methods are still better than those of the five other algorithms.For a low SNR = −5 dB, as shown in the bottom picture in figure 2, both the methods based on alternating optimization, namely ACDC and ACDC + LU , converge to local minima which are less desirable.The proposed algorithms are always able to converge to better results than the classical methods.The average numerical complexities and CPU time of all the algorithms over Monte Carlo trials are shown in table 2. It is observed that FFDIAG, LUJ1D and QRJ1D require a small amount of calculations, whereas ACDC + LU requires a large amount of calculations.The proposed JD + LU just costs a bit more flops and CPU time than ACDC, but it is still much more efficient.Concerning the JD + QR algorithm, it is more costly than JD + LU , with a comparable performance.We can then conclude that JD + LU offers the best performance/complexity compromise in these experiments.

Effect of SNR
In this section, we study the behaviors of the seven algorithms as a function of SNR.The dimensions of the 3-way array C N are set to N = 5 and K = 15.We repeat the experiments with SNR ranging from −30 dB to 50 dB with a step of 2 dB.The top picture in figure 3 depicts the average curves of α(A, Ã) of the seven algorithms as a function of SNR.The obtained results show that the performance of all the methods increases as SNR grows.For the unconstrained methods, generally ACDC performs better than FFDIAG, LUJ1D and QRJ1D.The nonnegativity constraint obviously helps ACDC + LU , JD + LU and JD + QR to improve the results for lower SNR values.The performance of ACDC and ACDC + LU remains stable for higher SNR values due to the small number of available sweeps and the lack of good initializations.Generally the proposed JD + LU and JD + QR algorithms outperform the others when SNR is between −20 dB and 30 dB and perform similar to FFDIAG, LUJ1D and QRJ1D when SNR is above 45 dB.The average numerical complexity and CPU time at each SNR level of all the methods in this experiment are shown in the bottom of figure 3. It shows that the proposed methods achieve better estimations of A and cost less flops and CPU time than ACDC + LU .The JD + LU gives the best performance/complexity trade-off for all the considered SNR values.

Effect of dimension K
In ICA, the third dimension K of the 3-way array C N ∈ R N ×N ×K corresponds to the number of covariance matrices at different lags, or the number of matrix slices derived from a cumulant array.In this section, we study the influence of K on the performance of the seven algorithms.The first and second dimensions of C N are set to N = 5.The SNR value is fixed to 10 dB.We repeat the experiment with K ranging from 3 to 55.The top picture in figure 4 shows the average curves of α(A, Ã) of all the algorithms as a function of K.For the five existing methods, ACDC, ACDC + LU , FFDIAG, LUJ1D and QRJ1D, their performance is quite stable with respect to K. The performance of the proposed methods progresses as K increases, and then practically stabilizes for high values of K.It indicates that after some point (e.g.K ≥ 20), the additional information brought by an increase of K does not further improve the results.The proposed JD + LU and JD + QR algorithms maintain competitive advantages through all the K values.The two images in the bottom of figure 4 present the average numerical complexity and CPU time of all the algorithms in this experiment, respectively.It shows that the numerical complexity of JD + LU and JD + QR is between that of ACDC and ACDC + LU .The JD + LU and JD + QR methods seem to be the most effective algorithms compared to the other methods.).Then the coherence ρ of D is defined as the maximum absolute cosine of angle ψ n,m between the columns of D as follows:

Effect of coherence of
The quantity ρ is also known as the modulus of uniqueness of JDC [43].By its definition (43), ρ falls in the range of [0, 1].The JDC problem is considered to be ill-conditioned when ρ is close to 1.Such an ill-conditioned problem can be met in ICA when A has nearly collinear column vectors.For example, in order to perform ICA, provided that all the sources are non Gaussian, which is often the case in practice, we can build a 3-way array C by stacking the matrix slices of the fourth order cumulant array of the observation data.Then the loading matrix D can be expressed as follows: where and where denotes the Khatri-Rao product.It can be observed that the coherence of the columns of A will influence the coherence of the matrix D. In the following test, the dimensions of the 3-way array C N are set to N = 5 and K = 15.The SNR value is fixed to 10 dB.In order to control ρ, firstly we randomly generate an orthogonal matrix D ∈ R 15×5 so that ρ = 0 by orthogonalizing a (15 × 5) random matrix.Secondly we rotate its 5 columns such that all the internal angles between any columns are equal to a predefined value ψ.Therefore ρ is only controlled by the angle ψ and equals to | cos(ψ)|.We repeat the experiment with the angle ψ ranging from 0 to π/2 with a step of π/60.A small ψ value means a large ρ value.The top picture in figure 5 displays the average curves of α(A, Ã) of all the algorithms as a function of ψ.It shows that the nonnegativity constrained methods ACDC + LU , JD + LU and JD + QR , outperform the unconstrained ones ACDC, FFDIAG, LUJ1D and QRJ1D.The proposed algorithms are more efficient, particularly when the coherence level is high.The average numerical complexity and CPU time displayed in the bottom of figure 5 indicate that the JD + LU algorithm provides the best performance/complexity compromise, while the JD + QR algorithm is also competitive with regard to ACDC + LU .

Effect of condition number of D (k)
When the JDC problem is considered, a diagonal matrix D (k) could contain some diagonal elements which, despite being non-zero, are many orders of magnitude lower than some other elements, leading to an ill-conditioned matrix C (k) .For the proposed methods, the inverse of such a matrix C (k) would contain numerical errors.
In this experiment, we study the performance of the seven algorithms as a function of the condition number of one of the diagonal matrices D (k) .The dimensions of the 3-way array C N are set to N = 5 and K = 15.The SNR value is set to 10 dB.We vary the condition number of the first diagonal matrix D (1) from 1 to 1000 by fixing the ratio of its largest diagonal element to its smallest diagonal element.The top picture in figure 6 displays the average curves of the estimating error α(A, Ã) of the seven algorithms as a function of the condition number of D (1) .The results reveal that a highly ill-conditioned diagonal matrix D (1) has a clear negative effect on the estimation accuracy of all the algorithms.The nonnegativity constrained methods ACDC + LU , JD + LU and JD + QR outperform the classical algorithms ACDC, FFDIAG, LUJ1D and QRJ1D whatever the condition number is.The proposed JD + LU and JD + QR algorithms maintain advantages when the condition number is less than 100.
Regarding the cases of larger condition numbers, ACDC + LU is more superior since it does not need to invert the highly ill-conditioned matrix.It is worthy pointing out that in practice, we can choose these sufficiently well-conditioned matrices C (k) for the proposed methods, whose condition numbers are below a predefined threshold.In addition, a weighted cost function for which weights would depend on the condition number of each matrix can be considered.On the other hand, the performance of the classical methods can also be improved by choosing a particular subset of available matrices [57] and by properly weighting the cost functions [49].In order to give a fair comparison, all the algorithms operate on the same set of matrices in all the experiments of this paper.In addition, the average numerical complexity and CPU time at each condition number of all the methods in this experiment are shown in the bottom of figure 6.It shows that the proposed methods give the best performance/complexity trade-off compared to ACDC + LU whatever the condition number is.
Test with a non-square matrix A As described in the section of practical issues, when a non-square matrix A ∈ R N ×P + with N > P is met in ICA, we propose to compress it by a nonnegative compression matrix W + ∈ R P ×N + [53], such that the resulting matrix Ā = W + A is a (P × P ) nonnegative square matrix.Then the proposed methods can be applied to estimate Ā. Similarly to the classical prewhitening, the nonnegative compression step could introduce numerical errors.In this experiment, we compare our methods to ACDC and ACDC + LU through a simulated ICA model.The latter algorithms can directly estimate a non-square matrix A from the fourth order cumulant matrix slices.The ICA model is established as follows: where zero-mean unit-variance source vector whose elements are independently drawn from a uniform distribution over T is the (N × 1) zero-mean unit-variance Gaussian noise vector, and A is the (N × 3) nonnegative mixing matrix whose components are independently drawn from a uniform distribution over [0, 1].In this context the SNR is defined by: For the proposed JD + LU and JD + QR algorithms, the given realization of {x[f ]} is compressed by means of a matrix W + ∈ R 3×N + computed using the method proposed in [53], leading to a 3-dimensional vector {x[f ]}.We compute the fourth order cumulant array of {x[f ]} and choose the first 3 matrix slices in order to build a 3-way array.Hence JD + LU and JD + QR decompose a (3 × 3 × 3) array.Once the compressed mixing matrix Ā is estimated, the original mixing matrix is obtained by equation (39).Regarding ACDC and ACDC + LU , the fourth order cumulant array of {x[f ]} is directly computed without compression.We apply ACDC and ACDC + LU on two 3way arrays with different third dimensions.The first array of dimension (N ×N ×3) is built by choosing the first 3 matrix slices from the fourth order cumulant array, while the second array of dimension (N × N × N ) is built using the first N matrix slices.We study the impact of the number of observations N on the performance of the JDC algorithms, by varying N from 4 to 24.The SNR value is fixed to 5 dB.The number of samples used to estimate the cumulants is set to 10 3 .Figure 7 shows the average curves of the estimating error α(A, Ã) of all the algorithms as a function of N .As it can be seen, when N ≤ 15, the larger the value of N , the more accurate estimation of A is achieved.When N > 15, the further increase of N does not bring significant improvement in terms of the estimation accuracy.ACDC and ACDC + LU give better results when the array with a larger third dimension is considered.Their results on (N × N × N ) arrays outperform the proposed methods when N = 4. ACDC + LU also gives the best estimation on (N × N × N ) arrays with N = 5.It suggests that the numerical errors introduced by the compression step limit the performance of the proposed methods when only a small number of observation is available.Such a negative effect can be partially compensated by using a large number of observations, since the proposed JD + LU and JD + QR methods maintain the highest performance in terms of estimation accuracy when N ≥ 6.The performance ACDC and ACDC + LU can be further improved by using a (N ×N ×N 2 ) array, which contains all the N 2 fourth order cumulant matrix slices.However, it leads to a higher numerical complexity especially for a large value of N .Regarding the proposed JD + LU and JD + QR methods, their performance can also be improved by using all the 9 matrix slices of the fourth order cumulant array of the compressed observation vector.Nevertheless, the experimental result has already shown that by using only a small number of matrix slices, JD + LU and JD + QR can maintain lower numerical complexities than ACDC and ACDC + LU , while achieving better estimating results, when a large value of N is considered.Therefore, despite the negative influence of the nonnegative compression, the proposed methods still offer a good performance/complexity compromise to estimate a non-square matrix A.

BSS application on MRS data
In this section, we aim to illustrate the potential capability of the proposed JD + LU and JD + QR algorithms for solving a real-life BSS problem by an application carried on simulated MRS data.
MRS is a powerful non-invasive analytical technique for analyzing the chemical content of 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 [58].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 (45), where x[f ] ∈ R N is the MRS observation vector, s[f ] ∈ 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 positive concentrations of the source metabolites.SNR is defined as in equation (46).In this experiment, two simulated MRS source metabolites {s 1 [f ]} and {s 2 [f ]}, namely the Choline (Cho) and Myo-inositol (Ins) (see figure 8 (b)), are generated by Lorentzian and Gaussian functions [59].Each of the sources contains 10 3 samples.The observation vector x[f ] is generated according to (45).The components of the (N × 2) mixing matrix A are randomly drawn from a uniform distribution.The additive noise ν[f ] is modeled as a zero-mean unit-variance Gaussian vector.The ICA methods based on the proposed JD + LU and JD + QR algorithms, namely JD + LU -ICA and JD + QR -ICA, consist of four steps: i) compressing {x[f ]} by means of a matrix W + ∈ R 2×N + [53], ii) estimating the fourth order cumulant array of the compressed observations and stacking all the cumulant matrix slices in a 3-way array, iii) decomposing the resulting 3-way array by means of JD + LU and JD + QR , respectively, and iv) reconstructing the sources.The JD + LU -ICA and JD + QR -ICA are compared to four state-of-the-art BSS algorithms, namely two efficient ICA methods CoM 2 [54] and SOBI [47], the NonNegative ICA (NNICA) method with a line search along the geodesic [55] and the NMF method [56] based on alternating nonnegativity least squares.The performance is assessed by means of the error α({s[f ]} T , {s[f ]} T ) between the true source s[f ] and its estimate s[f ], the numerical complexity, and the CPU time.To find out the detailed analysis of the numerical complexity of the classical ICA algorithms, the reader can refer to the book chapter [60].Figure 8 shows an example of the separation results of all the methods with N = 32 observations and a SNR of 10 dB.Regarding CoM 2 , SOBI, NNICA and NMF, there are some obvious disturbances presented in the estimated metabolites.As far as JD + LU -ICA and JD + QR -ICA are concerned, the estimated source metabolites are quasi-perfect.Furthermore, the comprehensive performance of all the methods will be studied by the following experiments with 200 independent Monte Carlo trials.
In the first experiment, the effect of the number of observations N is evaluated.The SNR is fixed to 10 dB.The six methods are compared with N ranging from 4 to 116 with a step of 4. The average curves of error α({s[f ]} T , {s[f ]} T ) as a function of N are shown in the left image of figure 9.It can be seen that the estimating errors of all the methods improve as N increases.It suggests that in noisy BSS contexts, using more sensors often yields better results.The proposed JD + LU -ICA and JD + QR -ICA methods maintain the competitive advantages.The average curves of the numerical complexities of this experiment are shown in the bottom left picture of figure 9. We can notice that the numerical complexities of all the methods increase with N .The complexities of JD + LU -ICA and JD + QR -ICA seem identical in the logarithmic scaled plot, that is because theoretically their complexities are mainly dominated by the computation of the nonnegative compression step and of the cumulants.Indeed, JD + LU -ICA is more computationally efficient than JD + QR -ICA in the step of CP decomposition of the cumulant array.This can be verified by the average CPU time of those methods, shown in the bottom right image of figure 9. We can observe that JD + LU -ICA is slower than CoM 2 , but it is faster than NNICA, SOBI and NMF.In the second experiment, we study the influence of SNR on the performance of the six methods.The number of observations N is set to 32.SNR is varied from 0 dB to 50 dB with a step of 2 dB.The average curves of the estimating error α({s[f ]} T , {s[f ]} T ), as well as those of the numerical complexities and CPU time as a function of SNR of all the six methods are shown in figure 10.The proposed JD + LU -ICA and JD + QR -ICA methods provide the best estimating results with moderate computational complexities and CPU time.Generally speaking, the JD + LU -ICA algorithm offers the best performance/complexity trade-off in this BSS experimental context.

Conclusion
We have proposed two methods, called JD + LU and JD + QR , in order to achieve the CP decomposition of semi-nonnegative semi-symmetric 3-way arrays.The nonnegativity constraint is imposed on the two symmetric modes of 3-way arrays by means of a square change of variable, giving rise to an unconstrained joint diagonalization by congruence problem.Therefore, the nonnegative loading matrix can be estimated by computing the joint diagonalizer.We consider the elementary LU and QR parameterizations of the Hadamard square root of the nonnegative joint diagonalizer, leading to two Jacobi-like optimization procedures.In each Jacobi-like iteration, the optimization is formulated into a minimization of a polynomial or rational function with respect to only one parameter.In addition, the numerical complexity for each algorithm has been analyzed.
The performance of the proposed JD + LU and JD + QR algorithms is evaluated with simulated semi-nonnegative semi-symmetric 3-way arrays.Four classical nonorthogonal JDC methods without nonnegativity constraints, including ACDC [23], FF-DIAG [25], LUJ1D [26], QRJ1D [26], and one nonnegative JDC method ACDC + LU [41], are tested as reference methods.The performance is assessed in terms of the matrix estimation accuracy, the numerical complexity, and the CPU time.The convergence property, the influence of SNR, the impact of dimension, the effect of coherence, and the influence of condition number, are extensively studied by Monte Carlo experiments.The obtained results show that the proposed algorithms offer better estimation accuracy by means of exploiting the nonnegativity a priori.The JD + LU algorithm provides the best performance/complexity compromise.The proposed algorithms are suitable tools for solving the ICA problems where a nonnegative mixing matrix is considered, such as in MRS.In this case, the 3-way array built by stacking the matrix slices of a HO cumulant array fulfills the seminonnegative semi-symmetric structure.We proposed two ICA methods, namely JD + LU -ICA and JD + QR -ICA, based on CP decomposition of the fourth order cumulant array using JD + LU and JD + QR , respectively.The source separation ability of the proposed algorithms is verified through a BSS application carried out on simulated MRS data.The JD + LU -ICA and JD + QR -ICA are compared to one NMF method [56], one nonnegative ICA method [55], and two classical ICA methods, namely CoM 2 [54] and SOBI [47].The performance is comprehensively studied as a function of the number of observations and of SNR.The experimental results 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.The JD + LU -ICA algorithm provides the best performance/complexity trade-off.
and an integer P , find a semi-nonnegative IND-SCAL model of C = [[A, A, D]], subject to the (N ×P ) matrix A having nonnegative components.
N ), (2, N − 1), . . .,(2,3), (1, N ), (1, N − 1), . . .,(1,2) Simulated semi-nonnegative INDSCAL model The synthetic semi-nonnegative semi-symmetric 3-way array C = [[A, A, D]] ∈ R N ×N ×K of rank N is generated randomly according to the semi-nonnegative IND-SCAL model (3).When used without further specification, all the algorithms are manipulated under the following conditions: i) Model generation: the loading matrix A ∈ R N ×N + is randomly drawn from a uniform distribution on the interval [0, 1].The loading matrix D ∈ R K×N is drawn from a Gaussian distribution with a mean of 1 and a deviation of 0.5.The pure array C is perturbed by a residual INDSCAL noise array V.The loading matrices of V are drawn from a zero-mean unit-variance Gaussian distribution.The resulting noisy 3-way array can be written as follows:

. The dimensions of the 3 -
way array C N ∈ R N ×N ×K are set to N = 5 and K = 15.The performance is assessed under three SNR conditions: SNR = −5, 10 and 25 dB, respectively.Figure 1 shows the convergence curves measured in terms of the cost function as a function of sweeps.It shows that FFDIAG, LUJ1D and QRJ1D exhibit fast convergence behavior.They converge in less than 20 sweeps.ACDC + LU decreases the cost function (4) quasi-linearly.ACDC and ACDC + LU do not converge in a maximum of 200 sweeps.The proposed JD + LU algorithm converges in about 100 sweeps when SNR = 25 dB and SNR = 10 dB, and it converges in about 40 sweeps when SNR = −5 dB.Regarding JD +

D
In this experiment, the effect of the coherence of the third loading matrix D of the 3-way array C = [[A, A, D]] is evaluated.Let d n and d m denote the n-th and m-th columns of D, respectively.The angle ψ n,m between d n and d m can be derived by using the following Euclidean dot product formula d T n d m = d n d m cos(ψ n,m

Figure 3
Figure 3 JDC performance versus SNR.The dimensions of C N are set to N = 5 and K = 15.Top: the average error α(A, Ã) evolution of all the algorithms as a function of SNR.Bottom: the average numerical complexities (left) and the CPU time (right) of all the algorithms, respectively.

Figure 4
Figure 4 JDC performance versus dimension K.The first and second dimensions of C N and the SNR value are set to N = 5 and SNR = 10 dB, respectively.Top: the average error α(A, Ã) evolution of all the algorithms as a function of dimension K. Bottom: the average numerical complexities (left) and the CPU time (right) of all the algorithms, respectively.

Figure 5
Figure 5 JDC performance versus coherence.The dimensions of C N and the SNR value are set to N = 5, K = 15, and SNR = 10 dB, respectively.Top: the average error α(A, Ã) evolution of all the algorithms as a function of internal angle ψ between any two columns of D. Bottom: the average numerical complexities (left) and the CPU time (right) of all the algorithms, respectively.

Figure 6
Figure 6 JDC performance versus condition number.The dimensions of C N and the SNR value are set to N = 5, K = 15, and SNR = 10 dB, respectively.Top: the average error α(A, Ã) evolution of all the algorithms as a function of the condition number of one of the diagonal matrices D(k) .Bottom: the average numerical complexities (left) and the CPU time (right) of all the algorithms, respectively.

2 s Table 2
Average numerical complexities (in flops) and computation time (in seconds) of the convergence experiment