Docteur De L'

Université De

Grenoble Alpes

Elaheh Sobhani

Thèse Dirigée Par

Pierre Comon

Christian Jutten

Massoud Babaie-Zadeh

M Laurent

Albera Professeur

M Eric

Moreau Professeur

M Reza Sameni

MARVASTI M Farrokh

Mme Sepideh Hajipour

M Pierre Comon

M Christian

M Massoud Babaie-Zadeh

Prof Pierre Comon

Prof Christian Jutten

Prof Massoud

Prof Farrokh Marvasti

Prof Laurent Albera

Prof Eric Moreau

THÈSE Pour obtenir le grade de

Keywords: Signal, Image, Parole, Télécommunication (SIPT) Tensor decomposition, Text mining, Unsupervised clustering, Hidden/latent variable, Third order moments, Forward-Backward Splitting, Proximal operator, Performance index, Permutation and scale ambiguity -----

Tensors or multi-way arrays are useful tools to identify unknown quantities thanks to the uniqueness of their decomposition. Tensor decompositions have been widely applied to obtain unknown components with physical meanings in many applications such as medical image and signal processing, hyperspectral images analysis, chemometrics, etc.

In this thesis, we investigate the application of tensor decomposition for probability estimations, which are required for some targeted data/text mining tasks such as unsupervised clustering of data/documents. Besides criticizing the existing tensor decomposition algorithms for probability estimations, we propose to apply some proper constrained tensor decompositions, which result in more reliable and accurate estimations. Moreover, we introduce an algorithm for constrained tensor decomposition, called Simple Forward-Backward Splitting (SFBS), which is based on Proximal Minimization. SFBS performs better than state-of-the-art in decomposing noisy tensors while computationally less expensive.

In addition, to evaluate the performance of tensor decomposition algorithms, we introduce an index that we name CorrIndex, which provides interpretable performance bounds, while keeping computational complexity to a reasonable level. Furthermore, we propose a method of moment estimation (standard averaging), which estimates the second and third order moments, with the same performance of state-of-the-art, but based on a much simpler concept, i.e. weighted averaging. Moreover, standard averaging performs better in small dimensions, and provides some advantages in terms of computational complexity.

3.2

The consistency of the moment estimates. The relative error in the estimation of second order moments (P ) and third order moments (T ) are plotted in solid and dotted line, respectively.

Note that the plot of standard averaging and Ruffini's estimator are superimposed. In this figure, the larger the size of the corpus, the fewer the number of corpuses, to keep the total block length (and hence the computational load) constant. . . List of Tables In each cell, top: non-negative AO-ADMM, bottom: Robust tensor power method. As it can be seen, the median and standard deviation of the error of non-negative AO-ADMM is always less than those of Robust tensor power method. . . . . 129

Comparison between all mentioned estimators in estimating

Performances in estimating probabilities φ and A in terms of

CorrIndex: the smaller the better. This experiment is carried out using a part of a well-known text data set, namely 20 Newsgroups. To be more precise, documents of the selected data set belong to four topics: "computer graphics", "baseball", "cryptography" and "Christianity". . The cumulative distribution of hidden variables:

List of Abbreviations and Acronyms

Φ(k) = Prob(h ≤ k)

Φ

The vector of cumulative distribution of topics:

[Φ(1), Φ(2), . . . , Φ(K) Ω A dictionary 1 N A vector of ones of dimension N A (n)
mode-n loading matrix (factor)

a k
The vector of conditional probabilities of dimension

D: [f k (1), f k (2), . . . , f k (D)] T A D×K The matrix of conditional probabilities: [a 1 , a 2 , . . . , a K ] A The

S P

The simplex set of dimension P t i 1 ,i 2 ,...,i n-1 ,:,i n+1 ,...,i N fiber mode-n of tensor T , where ":" denotes a MAT-LAB notation which means considering all values of way n Throughout this thesis, we consider a particular mixture model, called single-topic model [START_REF] David M Blei | Latent dirichlet allocation[END_REF], which can be used to describe a hidden relation among observed data (see Section 3.2 and Fig. 3.1 for more detailed definition of single-topic model). In a single-topic model, a set of multi-view (observed) variables are generated according to their corresponding hidden (latent or non-observed) variable. For instance, a corpus of documents can be described by a single-topic model: the topic of each document is seen as a hidden variable and the words of that particular document are its multi-view

T = λ; A (1) , A (2) , . . . , A (N ) CP decomposition of tensor T T (n) mode-n unfolding matrix of tensor T T = G; A (1) , A (2) , . . . , A (N )
1 CHAPTER 1. INTRODUCTION
variables generated according to its topic. Note that single-topic models are not limited to describe a corpus of documents, it can also be useful in modeling social networks [AGH + 14] and movie recommendation systems [START_REF] David M Blei | Latent dirichlet allocation[END_REF].

It has been shown that whenever observed data can be modeled by a single-topic model, some statistical parameters of the data are related to the tensor of moments [START_REF] Anandkumar | A method of moments for mixture models and hidden markov models[END_REF]. To be more precise, the probability of the hidden variable (topic) and the conditional probability of multi-view variables (words) can be estimated via the decomposition of the tensor of third order moments, which is a method of moments [START_REF] Pearson | Contributions to the mathematical theory of evolution[END_REF]. As mentioned before, by employing these estimated probabilities, the targeted data/text mining task (unsupervised clustering) can be performed by means of Bayes' law [START_REF] Ruffini | A new method of moments for latent variable models[END_REF].

The main goal of this thesis is to investigate the effect of the accuracy of tensor decompositions on estimated probabilities. In other words, we aim at answering the question of which type of tensor decomposition is more appropriate for probability estimation. In addition, exploring the drawbacks of the employed tensor decompositions for the above mentioned purpose in the literature and trying to propose proper solutions are other goals of this thesis.

The outcomes of the thesis are as follows:

(i) Regarding the estimation of the tensor of third order moments, we investigate the effect of moment estimators on the accuracy of probability estimations and the required size of corpus. Moreover, we propose an estimator which performs as well as the best estimator in the state-of-the-art but based on a much simpler concept, i.e. weighted averaging. In addition, our moment estimator performs slightly better for small dimension tensors and provides some advantages in terms of computational complexity.

(ii) Concerning tensor decomposition algorithms, we show experimentally that the probabilities estimated by algorithms employed in the literature are not always acceptable as probability values. More importantly, the employed tensor decomposition algorithms are not robust to additive noise and suffer dramatically from rounding errors, especially when working with a real corpus.

Nevertheless, we propose to consider proper constraints (such as nonnegativity or the simplex set) in tensor decompositions, which guarantees reliable results for the estimated probabilities, and provides much more accurate estimations. Moreover, these kinds of decompositions are more robust to additive noise and perform better on real text datasets.

In addition, we introduce an algorithm for constrained tensor decompositions based on the proximal concept, called Simple Forward-Backward Splitting (SFBS), which performs better than other constrained tensor decompositions especially in noisy scenarios and converges faster.

Further, we discuss how one can apply SFBS with sparsity constraints, which is needed for other applications such as co-clustering [START_REF] Evangelos | Coclustering as multilinear decomposition with sparse latent factors[END_REF].

Last but not least, due to utilizing a proximal regularization and the proximal concept, SFBS has a complete convergence guarantee even with non-convex constraints.

(iii) Due to permutation and scaling ambiguities in tensor decompositions, it is not straightforward to evaluate the performance of tensor decomposition algorithms. In this thesis, we introduce a performance index that we coin CorrIndex, which is invariant to permutation and scaling ambiguities. Contrary to existing methods in the literature, CorrIndex is more reliable in ill-conditioned cases, and also provides interpretable performance bounds, while keeping computational complexity to a reasonable level.

This thesis is organized as follows. In Chapter 2, we review some preliminaries about tensors and their decompositions.

Introduction

In this chapter, we bring the required notations and preliminaries in this report about the tensor and its decompositions. Moreover, the challenges and limitations of tensor decompositions in practice will be discussed along with the relevant methods to face these difficulties.

CHAPTER 2. TENSOR

Tensor: notations and preliminaries

Tensor Tensors can be considered as a multi-linear maps from a vector space to another one [START_REF] Comon | Tensors: a brief introduction[END_REF], or they are simply multi-way (multidimensional or multi-index) numerical arrays [CZPA09, CMDL + 15], or hypermatrices [START_REF] Lim | Tensors and hypermatrices[END_REF]. The order of an array refers to the number of its ways [START_REF] Comon | Tensor decompositions-state of the art and applications, keynote address in ima conf[END_REF]. Vectors and matrices are one-way and two-way arrays, respectively, but usually tensor refers to an array with three or more ways [SDLF + 17].

Products

• Tensor/outer product (⊗): It can be explained easily by considering an example. The tensor/outer product of a 3 rd order tensor, A, with entries A ijk and a 4 th order tensor, B, with entries B lmnp is a 7 th order tensor like

C = A ⊗ B whose components are C ijklmnp = A ijk B lmnp .
The orders add up under the tensor product.

• Contraction (• k ): Unlike the tensor product, which increases the order, the contraction decreases the sum of orders by 2 [START_REF] Comon | Tensors: a brief introduction[END_REF].

C = A • k B
indicates the contraction product over mode-k between two tensors, A and B, which is a summing up over the production of their mode-k entries. Note that if the order of A and B are D 1 and D 2 , C is a tensor of order D 1 + D 2 -2. For instance, in the previous examples of A and B, the contraction product over mode-2 is defined as follows:

C = A • 2 B → C iklnp = j A ijk B ljnp (2.1)
However, the contraction product between a matrix and a tensor of higher dimension is treated based on the order of putting them. To be more precise, G = A • k M is a contraction between the mode-k of A and the first matrix index of M , but

Q = P • k A is a contraction
between the mode-k of A and the second matrix index of P . This permits the compatibility with the usual matrix-matrix product.

CHAPTER 2. TENSOR

In the following, two examples of the contraction product over the mode-3 of a 3 rd order tensor, A, and two other matrices, namely M and P , can be found:

G = A • 3 M → G ijt = k A ijk M kt (2.2) Q = P • 3 A → Q ijs = k A ijk P sk (2.3)
• Kronecker product (⊠): Assume that A and B are two matrices of size I 1 × I 2 and I 3 × I 4 respectively. Their Kronecker product is a matrix C of size I 1 I 3 × I 2 I 4 defined by:

C = A ⊠ B def =        A 11 B A 12 B . . . A 1I 2 B A 21 B A 2I 2 B . . . A 2I 2 B . . . . . . . . . . . . A I 1 1 B A I 1 2 B . . . A I 1 I 2 B       
(2.4)

• Khatri-Rao product (⊙): This product is actually a column-wise Kronecker product, and it should be performed between two matrices with the same number of columns. If A I 1 ×I 2 = [a 1 , a 2 , . . . , a I 2 ] and

B I 3 ×I 2 = [b 1 , b 2 , . . . , b I 2 ]
, where a i and b i are the i th columns of the matrices A and B respectively, then:

C = A ⊙ B def = [a 1 ⊠ b 1 , a 2 ⊠ b 2 , . . . , a I 2 ⊠ b I 2 ] (2.5)
Unfolding or matricizing For convenience, tensors are sometimes transformed into matrices [START_REF] Comon | Tensors: a brief introduction[END_REF]. This transformation is called unfolding or flattering and can be done in each mode. The resulting matrix in mode n is called the mode-n unfolding [START_REF] Comon | Tensors: a brief introduction[END_REF] and is denoted by T (n) .

Unfolding can be executed in several manners, but in this thesis we shall use the one, which returns matrices with the same number of rows as in the chosen mode. The obtained vector by fixing all indices except one is called the fiber of tensor T I 1 ×I 2 ×.. The following example expresses the unfolding of the tensor A of dimension 2 × 2 × 2, depicted in Fig. 2.1, in its three modes:

A (1) =   A 111 A 121 A 112 A 122 A 211 A 221 A 212 A 222   A (2) =   A 111 A 211 A 112 A 212 A 121 A 221 A 122 A 222   A (3) =   A 111 A 211 A 121 A 221 A 112 A 212 A 122 A 222  

Exact tensor decomposition

In noiseless scenarios, an exact tensor decomposition is feasible. An advantage of exact tensor decomposition over exact matrix decomposition is the uniqueness of the decomposition under a mild condition (cf. Section 2.3.1).

Recall that the decomposition of a matrix as a summation of rank-1 matrices is not unique, for example if

M I×J = A I×K B T J×K = K k a k b T k , then M = (AQ)(Q -1 B T ) for any invertible matrix Q (QQ -1 = I).
CHAPTER 2. TENSOR

Canonical Polyadic (CP) tensor decomposition

A decomposable tensor of order N is a tensor product of N vectors [START_REF] Comon | Tensors: a brief introduction[END_REF],

i.e., D = a (1) ⊗ a (2) . . . ⊗ a (N ) . According to [START_REF] Frank | The expression of a tensor or a polyadic as a sum of products[END_REF], any tensor can be written as a linear combination of a finite number of decomposable tensors,

T = R r=1 λ r D(r), (2.6)
where T is a tensor of order N , and D(r) = a

(1)

r ⊗ a (2) r . . . ⊗ a (N )
r . The decomposition described in (2.6) is called the Polyadic decomposition [START_REF] Frank | The expression of a tensor or a polyadic as a sum of products[END_REF], and if it is unique, it is called the Canonical Polyadic (CP) decomposition or also CANDECOMP or PARAFAC [START_REF] Comon | Tensors: a brief introduction[END_REF]. In compact form, (2.6) can be represented as T = λ; A (1) , A (2) , . . . , A (N ) , in which λ is a coefficient vector of size R containing the values of λ r and a Rank For each tensor, the minimum value of R for which (2.6) holds is called the tensor rank [START_REF] Comon | Tensors: a brief introduction[END_REF]. Therefore, the rank of a decomposable tensor is one.

Uniqueness A sufficient condition of uniqueness of the CP decomposition is as follows [START_REF] Comon | Tenseurs en sciences des données[END_REF]:

2R + N -1 ≤ N n=1 krank{A (n) },
where N ≥ 3 and krank{.} denotes the Kruskal 1 rank of a matrix.

However, a necessary condition can be obtained according to the expected rank of a tensor for the uniqueness of its CP decomposition. Based on the CHAPTER 2. TENSOR number of knowns and unknowns in (2.6), the expected rank of a tensor is defined as [START_REF] Comon | Tensors: a brief introduction[END_REF]:

R • ≜ D 1 -N + N i=1 n i ,
where N is the order of the tensor of dimensions n 1 × . . . × n N and D = N i=1 n i . It has been shown in [START_REF] Chiantini | An algorithm for generic and low-rank specific identifiability of complex tensors[END_REF] that if D ≤ 15000 and R ≤ R • -1 ensures almost surely the uniqueness of decomposition (2.6). As mentioned before, if this polyadic decomposition is unique, it can be called the Canonical Polyadic (CP) decomposition [START_REF] Comon | Tensors: a brief introduction[END_REF].

Unfolding of CP Note that the mode-n unfolding of the tensor T can be expressed based on its CP decomposition, T = λ; A (1) , A (2) , . . . , A (N ) , as follows:

T (n) = A (n) diag(λ) A (N ) ⊙ . . . ⊙ A (n+1) ⊙ A (n-1) ⊙ . . . ⊙ A (1) T , (2.7)
where diag(λ) is a diagonal matrix by holding λ as its diagonal.

Tucker tensor decomposition

For any tensor of order N and of dimensions

n 1 × n 2 × . . . × n N , one can find N loading matrices A (1) , A (2) , . . . , A (N ) of dimension n 1 × R 1 , n 2 × R 2 , . . . , n N × R N respectively, and a core tensor, G, of dimension R 1 × R 2 ×
. . . × R N (smaller than or equal to dimension of T ) as follows:

T = A (1) • 1 A (2) • 2 . . . A (N ) • N G,
(2.8)

T = G; A (1) , A (2) , . . . , A (N ) .
(2.9)

If G in (2.8) is diagonal, then it expresses the CP decomposition, otherwise it is Tucker decomposition proposed by Tucker in 1966 [START_REF] Ledyard | Some mathematical notes on three-mode factor analysis[END_REF].

Uniqueness The uniqueness condition of the CP decomposition stated in Subsection 2.3.1 (R ≤ R • -1) can be refined by means of the dimensions

CHAPTER 2. TENSOR R 1 × R 2 × . . . × R N as follows [COV14, Com21]: R ≤ N i=1 R i 1 -N + N i=1 R i -1.
In addition, the rank R of the tensor T and the dimensions R 1 ×R 2 ×. . .×R N are related via the following inequalities [START_REF] Comon | Tenseurs en sciences des données[END_REF]:

max i {R i } ≤ R ≤ min i { j̸ =i R j }. Even if the minimal dimensions R 1 × R 2 × . . . × R N are used in (2.8),
Tucker decomposition is not unique. It can be seen easily by replacing each loading factor with its QR decomposition [START_REF] Comon | Tenseurs en sciences des données[END_REF], for example

A (i) = Q (i) R (i)
, where the matrices R (i) are of dimensions R i ×R i and

Q (i) H Q (i) = I.
Therefore, (2.9) can be rewritten as below:

T = G ′ ; Q (1) , Q (2) , . . . , Q (N ) , (2.10)
where the new core tensor G ′ has the following form: N ) .

G ′ = G; R (1) , R (2) , . . . , R (

2.3.3

Multi-Linear/Higher Order Singular Value Decomposition (MLSVD/HOSVD) Equation (2.10) can be seen as Tucker decomposition of T under the constraint that all loading factors should be semi-unitary (i.e.

Q (i) H Q (i) = I).
This decomposition is also known as Multi-Linear/Higher-Order Singular

Value Decomposition (MLSVD/HOSVD) [START_REF] Comon | Tensors: a brief introduction[END_REF].

Therefore, any tensor T admits a decomposition of MLSVD/HOSVD:

T = G; Q (1) , Q (2) , . . . , Q (N ) , (2.11)
where

Q (i) H Q (i) = I and G is a core tensor of dimension R 1 × R 2 × . . . × R N . Moreover, [R 1 , R 2 , . . . , R N ] is called the multi-linear rank of the tensor T [DLDMV00].

CHAPTER 2. TENSOR

Calculating exact MLSVD/HOSVD Assume that a multi-linear rank,

[R 1 , R 2 , . . . , R N ], has been chosen such that an exact MLSVD/HOSVD of T can be obtained. It can be shown that Q (i) in (2.11) is nothing else but the matrix of left singular vectors of the i th mode unfolding of T , i.e. T (i) :

T (i) = Q (i) Σ (i) V (i) H , (2.12)
and Σ (i) and V (i) are the diagonal matrix of singular values and the matrix of right singular vectors, respectively [START_REF] Comon | Tenseurs en sciences des données[END_REF]. In addition, G in (2.11) can be computed by the means of Q (i) H as follows:

G = T ; Q (1) H , Q (2) H , . . . , Q (N ) H .

Approximate tensor decomposition

Different types of rank (e.g. generic, typical) are discussed in [Com14, Section V.B]. The generic rank is the rank that is encountered with the probability one [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF][START_REF] Comon | Generic and typical ranks of multi-way arrays[END_REF], when the entries of a tensor are drawn independently according to a continuous probability distribution. In practice, almost always the tensor of data is corrupted by noise, which can be considered as an additive random tensor whose entries are random variables with continuous probability distributions. As a result, the rank of such a noisy tensor is generic [START_REF] Comon | Tensors: a brief introduction[END_REF]. This generic rank is much larger than the dimension, hence, its CP decomposition would not be unique.

Since the rank of the non-noisy part is much lower, it is typically preferred to perform a low-rank or a low multi-linear rank approximation, which not only eliminates partially the additive noise, but also results in less expensive computations. More importantly, the low-rank approximation may then have a unique CP decomposition.

Truncated MLSVD/HOSVD

Like truncated SVD in the matrix case, truncated MLSVD/HOSVD can be obtained by keeping a part of the columns of Q (i) corresponding to the largest singular values. To be more precise, considering truncated multi-linear rank

R t i ≤ R i , truncated MLSVD/HOSVD of T is its MLSVD/HOSVD by [R t 1 , R t 2 , . . . , R t N ].
In addition to being well-posed, truncated MLSVD/HOSVD is also useful as a pre-processing for other decompositions such as CP. These preprocessings include noise reduction (as it is a low multi-linear rank approximation) and dimension reduction (if the core tensor estimated via truncated MLSVD/HOSVD is used for the targeted decomposition).

Low-rank CP approximation

Most CP decomposition algorithms require the rank of tensor in advance.

Determining the rank in advance becomes more challenging in noisy scenarios, since the rank is generic. CORCONDIA is a well-known method for tensor rank estimation in terms of CP decomposition [START_REF] Bro | A new efficient method for determining the number of components in parafac models[END_REF]. In addition, as mentioned in Subsection 2.4.1, truncated MLSVD/HOSVD provides some pre-processings for targeted decompositions such as CP. Therefore, one can apply low-rank CP approximation on the core tensor estimated via truncated MLSVD/HOSVD, which has a lower rank compared to the primary tensor (due to the noise reduction) and has also lower dimensions, which result in less expensive computations. However, even if the rank of a tensor is estimated, the low-rank approximation is ill-posed [START_REF] Comon | Tensors: a brief introduction[END_REF]. In practice, in order to face this difficulty, some constraints must be added to CP decomposition, as described in the next Section.

Constrained CP decomposition

Since in practice, data stored in the form of a tensor are usually corrupted by noise, the best rank-R approximation must be estimated. Although lowrank approximation is useful and unavoidable, generally it is ill-posed [START_REF] Christopher | Most tensor problems are np-hard[END_REF][START_REF] De | Tensor rank and the illposedness of the best low-rank approximation problem[END_REF], since the set of tensors of rank at most R is not closed [START_REF] Comon | Tensors: a brief introduction[END_REF].

CHAPTER 2. TENSOR Therefore, imposing some constraints such as non-negativity [START_REF] Lim | Nonnegative approximations of nonnegative tensors[END_REF] or angular separation [START_REF] Lim | Blind multilinear identification[END_REF] in CP decomposition is proposed in the literature to overcome this difficulty.

We shall see the complete formulation of constrained CP decomposition in Section 5.5.1. Here, we bring a general description of some prevalent constraints in this context:

-Non-negativity: All the loading matrices and the coefficient vector are supposed to be non-negative even if the tensor to decompose is not completely non-negative due to noise. There is a rich literature on tensor decomposition under this constraint [CMDL + 15, JMC18].

-Simplex: This constraint usually appears in the models involving a probabilistic analysis [AGH + 14] (cf. Section 3.3.1.1). All or a part of the columns/rows of the loading matrices and/or the coefficient vector in (2.6) should belong to a simplex set defined by: S ≜ {x : x ≥ 0, ∥x∥ 1 = 1} -Orthogonality: it can be imposed between the columns of loading matrices or between decomposable tensors, D r , in (2.6) [START_REF] Comon | Tensors: a brief introduction[END_REF][START_REF] Kolda | Orthogonal tensor decompositions[END_REF].

The former constraint is widely used in blind source separation after the standardization stage [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF].

-Coherence: In order to guarantee the existence of a low-rank approximation, imposing orthogonality is sufficient but not necessary. In fact, imposing a constraint on the coherences of loading matrices is sufficient [SC15, [START_REF] Nazih | Using the proximal gradient and the accelerated proximal gradient as a canonical polyadic tensor decomposition algorithms in difficult situations[END_REF][START_REF] Lim | Blind multilinear identification[END_REF]. This constraint limits the coherence, i.e., the minimal angle between the columns of loading matrices [START_REF] Comon | Tensors: a brief introduction[END_REF], and this angle is generally much smaller than π/2. In addition, such an angular constraint often has a physical meaning, especially in antenna array processing [START_REF] Lim | Blind multilinear identification[END_REF].

-Sparsity: In some applications [START_REF] Hsieh | 2d sparse dictionary learning via tensor decomposition[END_REF], it may be needed to impose sparsity constraints on loading matrices or coefficient vectors. The CHAPTER 2. TENSOR exact sparsity constraint is also known as cardinality constraint, and defines a non-convex set [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF]. The cardinality of a vector can be measured by the ℓ 0 pseudo-norm (or the counting norm), which is the number of its non-zero entries. Since ℓ 0 is a non-convex function, sometimes its convex approximations such as the ℓ 1 norm [START_REF] Parikh | Proximal algorithms[END_REF], or Smoothed ℓ 0 (SL0) [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF] are used instead.

Conclusion

In this chapter, we reviewed some preliminaries about tensor and its decompositions in terms of exact and approximate decompositions. In addition, it was mentioned that imposing a proper constraint on the decomposition would be very efficient in some particular applications such as the one investigated in this thesis, i.e. data/text mining. In the next chapter, the relation between data/text mining with constrained tensor decomposition will be discussed in more details.

3 Data mining and tensors

Introduction

In this chapter, firstly, we review a particular data model, which can be used to describe a hidden relation among observed data. In addition, a task of data mining (clustering by means of estimated probabilities) based on this model will be described. Secondly, we shall explain an existing tensor approach as well as other points of view to perform this task. We will see the first step in utilizing the tensor approach is moment estimation, therefore, CHAPTER 3. DATA MINING AND TENSORS the last section of this chapter is devoted to discuss about different moment estimators and their properties. In this work, we assume the graphical model depicted in Fig. 3.1 as the multi-view model. In this model, a set of multi-view (observed) variables, {x 1 , x 2 , . . . , x L } are represented that are generated according to a conditional multi-nomial distribution conditioned to their corresponding hidden (latent or non-observed) variable h. This model is called the mixture of unigrams [START_REF] David M Blei | Latent dirichlet allocation[END_REF], since it can be obtained by augmenting unigram models with a random discrete hidden variable h, as the topic of each document.

Hidden

The model in Fig. 3.1 is also known as the single-topic model [RCG18, AGH + 14], since the multi-view variables (words in a corpus of documents)

are generated according to a single hidden variable (topic in each document).

However, there are also other more general models such as probabilistic Latent Semantic Indexing (pLSI) [START_REF] Hofmann | Probabilistic latent semantic indexing[END_REF] and Latent Dirichlet Allocation (LDA) [START_REF] David M Blei | Latent dirichlet allocation[END_REF], which attempt to correspond multi-view variables to several (not only one) hidden variables by means of a mixture of weights or Dirichlet distribution on a set of hidden variables. These models as well as other types of LVM such as Gaussian mixture and Hidden Markov Model (HMM) [START_REF] Michael | Graphical models[END_REF][START_REF] Anandkumar | A method of moments for mixture models and hidden markov models[END_REF] are out of the scope of this dissertation.

Employing a model as in Fig. 3.1 is very common in several tasks of data (e.g. text) mining [RCG18, AGH + 14, BNJ03]. In Section 3.3, these applications will be explained in brief.

Data mining with hidden and multiview variable models

The daily increasing amount of online available textual data have raised many challenging tasks such as unsupervised learning or clustering, which can be investigated by text mining [START_REF] Han | Data mining: concepts and techniques[END_REF]. Text mining as a special case of data mining is the center of attention of this thesis. Although the graphical model of Fig. 3.1 can be used for a wide range of learning problems [START_REF] Murphy | Machine Learning, A Probabilistic Perspective[END_REF],

in this section, we customize the formulation for the particular problem of text mining. The following explanation describes how one can encode textual data (words) into numerical vectors, hence applying matrix/tensor analysis on text become feasible. In addition, the required pre-processing steps to extract keywords from each text are briefly explained in Section 5.6.3, yet these steps are out of scope of this thesis.

Let L be the number of words in a given document, x ℓ the observed words, ℓ ∈ L = {1, 2, . . . , L}, and h a topic encoded into a discrete variable taking K possible integer values, say from H = {1, 2, . . . , K} with probability

φ(k) = Prob(h = k).
All words belong to a known encoded dictionary Ω = {u 1 , . . . , u D } of cardinality D. In other words, we can consider a mapping γ e (generally not injective) from L to {1, 2, . . . , D} such that x ℓ = u γe(ℓ) .

In the context of text mining, D would be the number of words and K the number of topics. The number of documents, N c , is generally larger than D.

Besides the probability of topics, the conditional probability of each word given the topic is also an important parameter in text mining tasks, e.g. unsupervised clustering of documents, as explained in the sequel. We denote the conditional probability of each word u d of the dictionary Ω, given a par-

ticular topic, h = k, by f k (d) = Prob(x = u d |h = k).
Therefore, according to the simplifying assumptions described in Section 3.3.1.1, the joint distribution of X = [x 1 , . . . , x L ] can be written as:

p X (u γe(1) , . . . , u γe(L) ) = K k=1 φ(k) f k (γ e (1)) . . . f k (γ e (L)), (3.1)
a k = [f k (1), f k (2), . . . , f k (D)] T Conditional probabilities of words to topic k A = [a 1 , a 2 , . . . , a K ]
Conditional probability of words which is referred to as the naive Bayes model. The main task is to estimate the quantities appearing in the right-hand side of (3.1) from realizations of X.

The term f k (γ e (ℓ)) in (3.1) represents the conditional dependence of the word x ℓ in a document to its topic k. In other words, the number of occurrences of a particular word in a document with a particular topic can be completely different from its occurrence in another document with another topic. For example, the word "win" is much more probable in a sport document than a document about cooking. On the other hand, the word "meat" may occur more in the latter document. For the convenience, above notations are listed in Table 3.1.

Estimating the probabilities in (3.1) is the first step in a variety of machine learning and data mining applications, some of which are described below.

CHAPTER 3. DATA MINING AND TENSORS

Unsupervised clustering [START_REF] Ruffini | A new method of moments for latent variable models[END_REF] In the following, as an example, it is explained how one can employ the estimated probabilities in order to cluster documents in an unsupervised way. However, this approach is applicable for other kinds of data admitting the graphical model in Fig. 3.1.

Assume that φ(k) and f k (γ e (ℓ)) in (3.1) for k = 1, . . . , K and ℓ = 1, . . . , L have been recovered. The hidden topic, h, of a document, say text, can be inferred according to the following conditional probability:

Prob(h = k|text) = Prob(text|h = k)φ(k) K k=1 Prob(text|h = k)φ(k) , (3.2)
and Prob(text|h = k), according to the simplifying assumptions described in Section 3.3.1.1, can be calculated as follows:

Prob(text|h = k) = Prob(u γe(1) , . . . , u γe(ℓ) , . . . , u γe(L) |h = k) = L ℓ=1 Prob(u γe(ℓ) |h = k) = L ℓ=1 f k (γ e (ℓ))
In fact, the document text is assigned to the topic k, which maximizes Prob(h = k|text). Note that the equality in (3.2) is Bayes' law.

Document classification [BNJ03]

The choice of feature is very important in some classification methods. In the document classification problem, although the individual words of all documents can be treated as a rich feature set, it is enormously large [START_REF] Joachims | Making large-scale svm learning practical[END_REF].

A solution would be a dimensionality reduction, which can be done by considering the estimated conditional probabilities in (3.1). In other words, f k (d) can be viewed as the d th feature of a document titled by topic k.

Moreover, since d refers to the index of a word in the assumed dictionary, the dimensionality reduction can be done, for instance, by keeping just two words of each document, which have the highest conditional probabilities.

After extracting features, document classification can be performed by applying a proper method such as Support Vector Machine (SVM) [START_REF] Steinwart | Support vector machines[END_REF].

CHAPTER 3. DATA MINING AND TENSORS

Movie recommendation [START_REF] David M Blei | Latent dirichlet allocation[END_REF] This application can be performed on a data set of users, which indicates their preferred movies. A portion of indicated movies are used to train the recommendation system, while a preferred movie of each user is held out for testing. The goal is to predict the held-out movie for each user, which can be viewed as a recommended movie in practice.

Considering users as the hidden variable and their preferred movies as the multi-view variables, one can assign the held-out movie by estimating the associated probabilities in (3.1), and then by following the procedure described in the above mentioned application of unsupervised clustering.

Tensor approach

In this section, the relation between the left-hand side of (3.1) and moments will be explained. In addition, two main tensor approaches for estimating the probabilities on the right-hand side of (3.1) will be described. As these methods utilize moments to estimate probabilities, they are considered as the method of moments [START_REF] Pearson | Contributions to the mathematical theory of evolution[END_REF].

Method of moments based on tensor decomposition

As mentioned before, the main goal is to estimate the probability of the hidden variable (topic) and the conditional probability of multi-view variables (words). In this section, we show how these probabilities can be estimated by means of moments and what is the role of tensor decomposition in this vein. To this end, some assumptions are required that are listed below [AGH + 14, AHK12a]:

• conditional probabilities do not depend on the order of words (exchangablility), e.g.

Prob(u γe(p) , u γe(q) , u γe(r) |h) = Prob(u γe(q) , u γe(r) , u γe(p) |h) CHAPTER 3. DATA MINING AND TENSORS

• words are conditionally independent given the topic, i.e.

(u γe(p) |h) ⊥ ⊥ (u γe(q) |h)

• words have the same conditional distribution given the topic, i.e.

(u γe(p) |h) ∼ (u γe(q) |h). Note that although these assumptions cannot be thoroughly satisfied in real textual data, such simplifying assumptions are required to identify probabilities through a tensor approach. In this thesis, we consider dictionaries and texts that satisfy these assumptions.

In the sequel, the second and third order moments will be needed:

P def = E x {x p ⊗ x q }, (3.3)
T def = E x {x p ⊗ x q ⊗ x r }, (3.4)
where P is a D ×D symmetric matrix and T a D ×D ×D symmetric tensor.

These moments do not depend on {p, q, r} provided these three integers are all different, which ensures the conditional independence assumed in (3.1).

But note that {γ e (p), γ e (q), γ e (r)} may not be different because γ e is not injective.

x ℓ is encoded to u d , and as in [AGH + 14], u d is chosen as the columns of the D×D identity matrix. For example, u d for the i th word in the dictionary is as follows:

u {d=i} (j) =    0 if j ̸ = i 1 if j = i , 1 ≤ j ≤ D.
Because of this choice for u d , these moments exhibit the following relations:

P = K k=1 φ k a k ⊗ a k , (3.5)
T = K k=1 φ k a k ⊗ a k ⊗ a k , (3.6)
where a k constructs the k th column of a matrix, say matrix A. Note that a k contains the values of f k (d) for all d and each k. We provide the proof of (3.5) here, which is not in the literature:

CHAPTER 3. DATA MINING AND TENSORS Proof.

P def = E {x p ⊗ x q } a = E h E x {x p |h} ⊗ E x {x q |h} b = K k=1 φ k D d=1 f k (d)u d ⊗ D d ′ =1 f k (d ′ )u d ′ c = K k=1 φ k a k ⊗ a k ,
where • a: the separation of E over (x, h) and using the conditional independency of x p , x q given h;

• b: the definition of E;

• c: the specific encoding of u d into the columns of the D × D identity matrix.

The proof of (3.6) is exactly similar.

Rewriting (3.1) as (3.3) and (3.4) reveals the nice property that arrays P and T are actually the joint probabilities of observations, i.e., P ij =

Prob{x p = u i , x q = u j } and T ijk = Prob{x p = u i , x q = u j , x r = u k }.
Once the empirical approximation of the moments in (3.3) and (3.4) is obtained, the probability of the hidden variable (topic) and the conditional probabilities of multi-view variables (words), i.e. φ and A, can be estimated via tensor decomposition according to (3.6). As we shall see in Chapter 5, this decomposition may be either developed in a constrained manner (under non-negative constraints as described in Section 5.5.2.1 or under simplex set constraints as marked out in Section 5.5.2.2) or be handled in an unconstrained manner as will be expressed in Section 5.4. In addition, it may be performed with or without utilizing the second order moments matrix, P .

The advantages and disadvantages of each method are investigated experimentally in Section 5.6.2.

CHAPTER 3. DATA MINING AND TENSORS

Diagonalization of two moment matrices

In [START_REF] Anandkumar | A method of moments for mixture models and hidden Markov models[END_REF], a joint diagonalization algorithm is promoted uses two moment matrices, namely P and a matrix obtained by the contraction T (η) = T • 3 η, where η is a randomly drawn vector. The idea may seem interesting, but the algorithm unfortunately has never been implemented and tested, according to the available literature.

The Algorithm A of [START_REF] Anandkumar | A method of moments for mixture models and hidden Markov models[END_REF] starts with an SVD of P as P = U ΣV T .

Then, the matrix B(η) = U T T (η)V (U T P V ) -1 is computed, as well as its K dominant eigenvectors, ξ k . Then, an estimate of the columns of matrix A (the conditional probabilities of multi-view variables) is given by a k = U ξ k , up to a scaling factor depending on how a k are normalized. The "eigenvalues" of T can be obtained in a second stage by contraction as

φ k = T • 1 a k • 2 a k • 3 a k .

Other methods and points-of-view

Expectation Maximization (EM) For decades, the most popular unsupervised and heuristic approach to learn the parameters of a LVM model was Expectation Maximization (EM) [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF], which is an iterative local search based on Maximum Likelihood (ML) [AGH + 14].

Although EM has several merits such as the ease of understanding, implementing and utility for any LVM, it suffers from slow convergence specially when the model dimension increases and also from suboptimal local results [START_REF] Ruffini | A new method of moments for latent variable models[END_REF]. Moreover, not only it requires users to apply many heuristics to obtain acceptable results, but also for some kinds of models such as Latent trees, it is NP-hard [AGH + 14].

In addition, supposing a prior (e.g. Dirichlet) is needed in some particular LVM -such as the LDA model [START_REF] David M Blei | Latent dirichlet allocation[END_REF] -and some existing EM methods [START_REF] Nigam | Text classification from labeled and unlabeled documents using em[END_REF], and increasing the generality of the model is more costly and challenging in terms of unsupervised estimation procedure [START_REF] Anandkumar | A method of moments for mixture models and hidden Markov models[END_REF].

Latent Semantic Analysis (LSA) [DDF + 90] Analysis of texts and documents is a great challenge in Natural Language Processing (NLP), Information Retrieval (IR) and information filtering [START_REF] Hofmann | Probabilistic latent semantic indexing[END_REF].

One of the problems in organizing, searching and understanding such a vast amount of digitalized documents is learning the meaning and the usage of the words with a data-driven model [START_REF] David | Probabilistic topic models[END_REF]. In other words, one of the big challenges is the transition from the lexical level (i.e. actual written text) to the semantical level (i.e. the intention of the writer) [START_REF] Hofmann | Probabilistic latent semantic indexing[END_REF]. Some issues such as polysems (i.e. a word with several meanings) and synonyms (i.e. different words with similar meanings) make this challenge even harder.

Latent Semantic Analysis (LSA) [DDF + 90] is a well-known method, which maps the high dimensional space of the words of the vocabulary into a lower dimensional space, called latent semantic space. This method works with a matrix, W , whose rows and columns are referred to the documents of the corpus and the words of the vocabulary, respectively. W (i, j) shows how often the j th word occurs in the i th document, and the inner product of two columns of W shows how much those two words are co-occurred 2 For instance, "red" is hyponym of "color", and "color" is hypernym of "red", and "red" is co-hyponym of "blue".

CHAPTER 3. DATA MINING AND TENSORS types of latent relation.

Unlike the method of moments, which tries to estimate some latent statistics and the parameters of the observed data, LSA and its variants try to measure semantic data relation such as synonym, antonym, hypernym, etc.

In addition, LSA and its variants are based on cosine similarity of the raw observed data, whereas the method of moments attempts to estimate parameters using the empirical moments of the observed data. In general, LSA and its variants are applicable to understand the semantical layer from the lexical level, while the method of moments are useful to extract the informative features of each of documents in the corpus.

The approach used in this thesis In this thesis, we focus on the twophase tensor decomposition approach for data mining (i.e. estimating the probabilities) with the graphical model in Fig. 3.1, which was described in section 3.3.1.1. In order to make the procedure more clear, we bring a brief review below:

In the first phase, the second and third order moments are approximated by observing the multi-view variables (words). Then, in the second phase, one tries to estimate the probabilities by decomposing or using the tensor of third order sample moments (sometime the second order sample moments is also used).

The more precise the sample moments are approximated in the first phase, the more accurate estimation of the probabilities is obtained in the second phase. We shall see the effect of moments approximation in Section 5.6.2 in practice. In addition, in Section 3.4, some methods of moments approximation are reviewed.

Moment estimation

As mentioned before, the first step in the procedure of data mining (the probability estimation) by means of tensor decomposition described in Sec-CHAPTER 3. DATA MINING AND TENSORS tion 3.3.1.1 is the moment estimation. In this section, first, two generative processes are expressed, then, some of the methods of the moment estimation are reviewed (Section 3.4.3) and are introduced (Section 3.4.4). Finally, these moment estimators are compared by some simulations in Section 3.4.5.

It is worth first reviewing some notations of Section 3.3:

• L: the number of words in a given document,

• x ℓ , ℓ ∈ L = {1, 2, . . . , L}: the observed words in a document,

• h: a topic encoded into a discrete variable taking K possible integer values, say in H = {1, 2, . . . , k, . . . , K} with the probability φ(k) =

Prob(h = k) • Ω = {u 1 , . . . , u d , . . . , u D }: the encoded dictionary of cardinality D • f k (d) = Prob(x = u d |h = k): the conditional

Generative process

In this subsection, we explain how one can generate a set of synthetic encoded words that are related to the same topic (cf. Fig. 3.1), with these two properties simultaneously:

• conditionally independent given the topic

• having the same conditional distribution

In the following, we describe two generative processes, which generate a topic (as hidden variable) and words (as multi-view variables) such that the above mentioned properties will be developed. • draw z ∈ [0, 1] with uniform distribution, and pick a topic (hidden variable), h = Φ -1 (z), by selecting the first entry in Φ that is larger than z.

• assuming that h = k has been drawn in the first step, then, for each ℓ ∈ {1, 2, . . . , L},

draw z ℓ ∈ [0, 1] with uniform distribution, and pick a word (multiview variable), d = γ e (ℓ) = F -1 k (z ℓ ), by selecting the first entry in the k th column of F that is larger than z ℓ .

set x ℓ = u γe(ℓ) = u d .

Generative process based on the multinomial distribution

This generative process is stated in [START_REF] Ruffini | A new method of moments for latent variable models[END_REF]. After defining the distributions φ(k) and f k (d) for all (k, u d ) ∈ H × Ω, a topic and words are drawn based on some multinomial distributions with event probabilities φ and A, respectively. The "multinomial distribution" was firstly introduced by Ronald Fisher in 1925 [START_REF] Fisher | Theory of statistical estimation[END_REF]. Generally speaking, it is an extension of the binomial distribution to an event with more than two outcomes [START_REF] Upton | A Dictionary of Statistics[END_REF]. Some required moments of the multinomial distribution are expressed in Remark 1.

The details of this generative algorithm are as follows:

CHAPTER 3. DATA MINING AND TENSORS

• draw the topic (hidden variable), h, according to a multinomial distribution3 with K event probabilities of φ K×1 .

• assuming that h = k has been drawn in the first step, then, for each ℓ ∈ {1, 2, . . . , L}, -draw a word (multi-view variable), d = γ e (ℓ), according to a multinomial distribution with D event probabilities of a D×1 k .

set x ℓ = u γe(ℓ) = u d .

Moment consistency

Reminding that second and third order moments in (3.3) and (3.4) are defined as follows:

P def = E {x p ⊗ x q }, T def = E {x p ⊗ x q ⊗ x r }.
As it is shown in Section 3.3.1.1, by encoding x ℓ to u d and by choosing u d as the columns of the D × D identity matrix, the above mentioned moments have the following matrix and tensor structure (mentioned in (3.5) and (3.6)):

P = K k=1 φ k a k ⊗ a k , T = K k=1 φ k a k ⊗ a k ⊗ a k .
These moments, using the distributions φ k and f k (d), can be considered as true moments, while any sample estimate of (3.3) and (3.4) should converge to the true moments P and T defined in (3.5) and (3.6). Hence, it is essential to evaluate the consistency of the moment estimator in order to make sure that sample moments indeed converge to the true joint probabilities. It is preferred that this convergence occurs by observing the least possible number of realizations of x ℓ .

CHAPTER 3. DATA MINING AND TENSORS

To this end, first, sample moments are calculated via a moment estimator, then, the difference between true moments ((3.5) and (3.6)) and the estimated ones will be reported in terms of Euclidean norm. Therefore, the moment consistency can be considered as a measure of evaluating the performance of moment estimators, and we shall compare some algorithms of the moment estimation in Section 3.4.5 in terms of moment consistency.

State-of-the-art

Simple averaging [AFH + 12]

This estimator is simply the averages of large number of realizations (a corpus of documents), say N c , to obtain an acceptable approximation of P and T . In each realization, a random encoded topic, k, is drawn in the way described in the first step of generative processes in Sections 3.4.1.1 and 3.4.1.2. Then, according to the chosen topic, three random encoded words (it is assumed that a document consists of at least three words) are drawn, {x p = u γe(p) , x q = u γe(q) , x r = u γe(r) }, based on the second step of generative processes mentioned in Sections 3.4.1.1 and 3.4.1.2. At the end, by using the following averages, an empirical approximation of the second and the third order moments are obtained with sample moments below:

P e = 1 N c Nc n=1 u γe(p) ⊗ u γe(q) , (3.7
)

T e = 1 N c Nc n=1 u γe(p) ⊗ u γe(q) ⊗ u γe(r) . (3.8)
Note that in this method of estimation, the length of each document (L n )

does not matter, and just three words (can be the first, the middle and the last word) in each document participate in the sample moment estimation.

Due to these facts, providing a reasonable moment consistency with simple averaging is very difficult and requires a large size corpus, i.e. large N c in (3.7) and (3.8). We shall show this drawback of simple averaging in practice by some simulations in Section 3.4.5.

CHAPTER 3. DATA MINING AND TENSORS

Zou's estimator [ZHPA13]

For the rest of this section, a well-known concept in text mining, called bagof-words, is required, which is a representation of each document. Assume that a document consists of L n words (x ℓ , ℓ = 1, . . . , L n ) that are encoded to u d = u γe(ℓ) . Then, the bag-of-words of such a document is a vector of dimension D, and is defined as follows:

b n def = Ln ℓ=1 u γe(ℓ) , (3.9)
where b n (d) shows how many times the d th word appears in the document n.

The Zou's estimator (P D×D Zou , T D×D×D

Zou

) is as follows:

P Zou (i, i) = 1 N c Nc n=1 b n (i) 2 -b n (i) L n (L n -1) (3.10) P Zou (i, j) = 1 N c Nc n=1 b n (i) -b n (j) L n (L n -1) (3.11) T Zou (i, i, i) = 1 N c Nc n=1 b n (i) 3 -3b n (i) 2 + 2b n (i) L n (L n -1)(L n -2) (3.12) T Zou (i, i, j) = 1 N c Nc n=1 b n (i) 2 b n (j) -b n (i)b n (j) L n (L n -1)(L n -2) (3.13) T Zou (i, j, k) = 1 N c Nc n=1 b n (i)b n (j)b n (k) L n (L n -1)(L n -2) , (3.14)
where N c is the number of documents in the considered corpus. Equations

Ruffini's estimator [RCG18]

With the same assumptions of Zou's estimator, Ruffini's estimators, i.e. P Ruffini and T Ruffini , are defined as follows (i < j):

P Ruffini (i, j) = Nc n=1 b n (i)(b n (j) -δ i=j ) Nc n=1 L n (L n -1) (3.15) T Ruffini (i, j, k) = δ i<j<k Nc n=1 b n (i)b n (j)b n (k) Nc n=1 L n (L n -1)(L n -2) + δ (i=j<k)∨(i<j=k) Nc n=1 b n (i)(b n (j) -1)b n (k) Nc n=1 L n (L n -1)(L n -2) + δ i=j=k Nc n=1 b n (i)(b n (j) -1)(b n (k) -2) Nc n=1 L n (L n -1)(L n -2) . (3.16)
By employing bag of words of each document in the corpus, nominators of the above expressions correspond to implementing moments (summation and production of random variables), and their denominators are the coefficients, which make this estimator unbiased.

Note that by computing the upper triangular part of P Ruffini and T Ruffini , and then, by setting the lower triangular part identical to the upper side, symmetrization would be done.

Like Zou's estimator, all the words of the documents in a corpus along with the length of each document have been also considered in Ruffini's estimator. However, in Zou's estimator, the averaging is done for each document, and then, all the documents are averaged together with the same weight. But Ruffini's estimator averages all the documents according to a particular weight for each of them, i.e.

L n (L n -1) Nc n=1 L n (L n -1) [RCG18]
. Therefor, if all the documents have the same length (L n = L), the two estimators will produce the same estimation. Moreover, Ruffini's estimator is less sensitive to the noise [START_REF] Ruffini | A new method of moments for latent variable models[END_REF].

It is shown that Ruffini's estimator is unbiased [START_REF] Ruffini | A new method of moments for latent variable models[END_REF]. We make this CHAPTER 3. DATA MINING AND TENSORS proof easier to follow by providing much more details below (i ̸ = j):

E{P Ruffini (i, j)} = E Nc n=1 b n (i)b n (j) Nc n=1 L n (L n -1) = 1 Nc n=1 L n (L n -1) Nc n=1 E{b n (i)b n (j)} = 1 Nc n=1 L n (L n -1) Nc n=1 K k=1 φ(k)E{b n (i)b n (j)|h = k} = 1 Nc n=1 L n (L n -1) Nc n=1 K k=1 φ(k)[Cov(b n (i)b n (j)|h = k) I + E(b n (i)|h = k)E(b n (j)|h = k) II ] ( * ) = 1 Nc n=1 L n (L n -1) Nc n=1 K k=1 φ(k)[-L n a k (i)a k (j) I + L n a k (i)L n a k (j) II ] = K k=1 φ(k)a k (i)a k (j) = P (i, j),
where "Cov(., .)" denotes the covariance of two random variables, and the equality ( * ) is due to the properties of multinomial distribution, which are described in Remark 1.

Remark 1. Suppose that y is a random vector of dimension D holding multinomial distribution, i.e. y ∼ multinomial (t, [p 1 , . . . , p D ]), where "t" is the number of the trials of "D" events with probabilities [p 1 , . . . , p D ], and of course D d=1 p d = 1. According to the multinomial distribution [UC08], we have: Save the index of words of document "n".

• mean = E{y(d)} = tp d , • variance = E{y(d), y(d)} = tp d (1 -p d ), • covariance = E{y(d), y(d ′ )} = -tp d p d ′
4:

P temp = zeros(D,D) 5:
count-pair = 0

6:
for all possible pairs of words, ex. (word 1 , word 2 ) do 7:

P temp (word 1 , word 2 ) = 1 + P temp (word 1 , word 2 ) 8: count-pair = 1 + count-pair 9:
end for 10:

P Stand-Ave = P Stand-Ave + L n Nc n=1 L n P temp count-pair 11: end for 12: Symmetrize P Stand-Ave

Proposed: Standard averaging

By taking into account the length and all the words (not only three) of each document, we propose an algorithm, which we coin as standard averaging, and can be seen as a weighted averaging based on the length of each document.

Algorithm 1 describes standard averaging in details. Note that symmetrization in line 12 is executed by computing P + P T 2 . We observed experimentally that this method of symmetrization shows a better performance compared to that of Ruffini's estimator in which the values of the lower triangular part will be set identical to the upper one at the end of the algorithm. The algorithm for the third order moments is similar, except that in line 6, all possible triples would be considered. Note that the symmetrization of a third order tensor (cf. line 12) requires averaging over six permutations of indices.

CHAPTER 3. DATA MINING AND TENSORS

Simulations in Section 3.4.5 show that standard averaging performs better than the most recent estimator, i.e. Ruffini's estimator, when the dimension of dictionary, D, is not very large. More details about performances will be given in Section 3.4.5. Moreover, standard averaging provides some advantages in terms of computational complexity:

• standard averaging requires less number of multiplications in the estimation of T ,

• the number of required additions in estimating P does not depend on the size of dictionary (D), which is an advantage in case of large dimension dictionaries.

Comparison of computational complexity in terms of number of multiplications and additions are explained in the next paragraph.

Computational complexity By assuming the same method of symmetrization, we can compare the computational complexity of standard averaging and Ruffini's estimator in terms of multiplications and additions. In this calculation, we consider the symmetrization of Ruffini's estimator, i.e. setting the lower triangular part identical to the upper one. Note that to symmetrize the output of Algorithm 1 in this way, the symmetric pair of (word 1 , word 2 ),

i.e. (word 2 , word 1 ), will be ignored in the set of all possible pairs (line 6 of Algorithm 1). Remind that, as P and T are symmetric, the number of their free parameters are D(D + 1) 2 and D+3-1

3

, respectively [START_REF] Comon | Tensors: a brief introduction[END_REF].

One can observe that the number of required multiplications for comput-

P T Ruffini (1 + N c ) D(D + 1) 2 + N c (1 + 2N c ) D+2 3 + 2N c standard averaging N c D(D + 1) 2 + 2N c N c D+2 3 + 2N c 2N c +2 Nc n=1 Ln 2 .
Therefore, the number of additions in standard averaging does not depend on D, which could be preferable for large values of D. These results are also mentioned in Table 3.2.

Simulations

All computer experiments reported in this section have been executed either on a laptop with a processor of 3.1 GHz Intel Core i5, 16 GB RAM or on a PC with a processor of 3.2 GHz Intel Core i5, 8 GB RAM, both running macOS Mojav and MATLAB 2019a.

The effect of corpus size on the moment consistency

As mentioned before, simple averaging does not provide acceptable moment consistency and a naive way to compensate this drawback is increasing the size of the corpus (the number of documents or realizations of the topic/hidden variable and words/multi-view variables).

In order to show this challenge in practice, we generated N c realizations of our synthetic data set, x ℓ = u γe(ℓ) , (according to the method explained in of Euclidean norm. As can be seen in Fig. 3.2, simple averaging estimation, in the way described in Section 3.4.3.1, is consistent, if it performs on a corpus of size N c > 10000. Ruffini's estimator and standard averaging have the same relative error (they are superimposed in Fig. 3.2), however, we shall see the difference between the performance of Ruffini's estimator and of standard averaging in Section 3.4.5.2 by investigating the effect of dictionary cardinality.

The effect of dictionary cardinality on the moment consistency

In this section, we investigate the effect of dictionary cardinality, D, on the performances of the moment estimators. For all the experiments, we con- (the number of words) or K (the number of topics). We shall conclude that standard averaging performs better than Ruffini's estimator, when dictionary cardinality, D, is less than 50. As it can be seen in Figs. 3.3 and 3.4, not only standard averaging is a much more better estimator than simple averaging, but also it performs slightly better than Ruffini's estimator when dictionary cardinality, D, is less than 50.

In order to show that the only critical parameter is dictionary cardinality and not K, we performed two other experiments in Fig. 3

negligible).

As mentioned before, the length of a document, L, is a critical factor in improving the performance. When dictionary cardinality is small, the influence of document length is more significant, since the related words depending on the topic are much more repetitive in such a document. Ruffini's estimator decreases this influence by the denominator in its formula, e.g. by dividing by Nc n=1 L n (L n -1), while standard averaging is simply a weighted averaging with respect to the length of all documents (e.g. by dividing over Nc n=1 L n ). Therefore, standard averaging is able to better show the effect of the length, L n , in case of small D.

Conclusion

In this chapter, we described the hidden and multi-view variable model, which can be fitted to several multi-modal events. Then, we explained a task of data mining, i.e. probability estimation, which can be followed to perform an ultimate data mining task such as an unsupervised clustering.

In Section 3.3, we described a tensor approach for this probability estimation as well as some other existing methods and points of view. It is mentioned that the first step in utilizing the tensor decomposition to approximate the probabilities of hidden variables and the conditional probabilities of multi-view variables is moment estimation. In Sections 3.4.3 and 3.4.4, we described some existing methods of moment estimation as well as our proposed method, standard averaging, along with some simulation comparisons in Section 3.4.5, using a synthetic generative model.

The next two chapters are devoted to the second step of the tensor approach for data/text mining, i.e. the constrained tensor decomposition and the performance index needed to evaluate different decomposition algorithms.

Introduction

As mentioned in Section 3.3.1, tensor decomposition is the second step in estimating probabilities from observed data. In order to compare different methods of decomposition, a proper performance index is needed to be applied on the estimated coefficient vector and loading matrices. Due to permutation and scaling ambiguities, evaluating these variables is not straightforward. Section 4.2 explains the existing challenges with this issue.

In Section 4.3, the problem is formulated and a review of previous methods and measures is provided. The proposed index, CorrIndex, is presented in

Challenges in measuring performance

Permutation and scaling ambiguities are relevant issues in applications such as tensor decomposition [START_REF] Comon | Tensors: a brief introduction[END_REF] and Blind Source Separation (BSS) [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF].

Scaling refers to multiplication by a diagonal matrix with non-zero entries, which may be complex in the most general case, and permutation refers to the column permutation of a matrix, which is equivalent to multiplication by a permutation matrix. Firstly, these two ambiguities are inherent in tensor representations, by definition of tensors [START_REF] Comon | Tensors: a brief introduction[END_REF]. Secondly, in BSS, statistical independence is not affected by scaling or permutation of the sources [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]. A mixing (or demixing) matrix can then only be estimated up to these ambiguities under the independence assumption. Although it is impossible to eliminate these ambiguities when working with real data sets, where the original parameters are not available, it is feasible to overcome these uncertainties in evaluating algorithm performance on synthetic data sets. Furthermore, reasonable comparisons on synthetic data sets are very helpful to choose adequately an appropriate algorithm to be applied on real data sets. Therefore, in order to report reasonably the performance indices of existing algorithms on synthetic data sets where the desired parameters are accessible, it is important to employ proper methods to measure the performances.

Assume that the original and estimated components have been normalized, then the only remaining ambiguities are the permutation and scaling with complex numbers of unit modulus. The existing approaches to measure the performances of the algorithms of BSS and tensor decomposition can be classified in three main categories: "greedy approaches", "graph-based methods" and "invariant indices". Greedy approaches [FIW + 20, CLDA09, CKAC14, BBK18] try to assign the most correlated components estimated by an algorithm, and then compute the error of estimation or decomposition.

Although most of these methods return back an estimated permutation as well as a performance index, they are not reliable in noisy conditions. In other words, the reported index (which is a criterion to evaluate the performance of algorithms in estimating components) by these kinds of methods depends directly on the manner of computing and analyzing the correlation matrix.

Graph-based methods [START_REF] Harold W Kuhn | The hungarian method for the assignment problem[END_REF][START_REF] Galil | Efficient algorithms for finding maximal matching in graphs[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF] are originated from the well-known optimal assignment problem [START_REF] Duan | Linear-time approximation for maximum weight matching[END_REF], which is itself a particular case of the optimal transport problem [PC + 19]. Although these kinds of methods have the guarantee to find the optimal permutation, they are computationally expensive (as we shall see, the minimum cost is approximately 8N 3 flops with a correlation matrix of size N × N ), especially when the correlation matrix is large.

However, the viewpoint of a third category, namely invariant indices [Com94, MM94, ACY + 96], differs from the latter approaches. These invariant indices measure the performance regardless of permutation and scaling, and yield an index that can directly be used to compare algorithms. The reported indices of [Com94, MM94, ACY + 96] are invariant to permutation and scaling, and the index of [START_REF] Comon | Independent component analysis, a new concept? Signal Proc[END_REF] provides the guarantee of a zero distance between estimated and original matrices up to column permutation and scaling, when the obtained index is zero. Nevertheless, the index of [START_REF] Comon | Independent component analysis, a new concept? Signal Proc[END_REF] is not bounded from above. More importantly, the upper bounds of indices of [MM94, ACY + 96] have not been investigated, and it seems that these bounds are not easy to interpret. In addition, these methods are in the literature of source separation, and the indices introduced therein utilize the inverse (or pseudo-inverse) of the mixing matrix, which may involve an additional computational burden.

In next section, these three kinds of methods for measuring the perfor-CHAPTER 4. PERFORMANCE INDEX mance along with the problem formulation will be explained. Λ is a diagonal matrix with unit modulus entries and W is the error of estimation A, which can be modeled as an additive noise. Since relevant issues in considered applications are permutation and scale ambiguity of columns of A comparing to A, then right-multiplication of P σ Λ is utilized here. More formally, the goal is to measure the gap defined below:

State-of-the-art

Let A = [a 1 , a 2 , . . . , a N ] ∈ C M ×N and A = [â 1 , â2 , . . . , âN ] ∈ C M ×N
ϵ 0 (A, A) = min σ,Λ ∥AP σ Λ -A∥ 2 F (4.1)
This gap can be computed with or without estimating permutation σ explicitly. Seeking the optimal permutation σ can be written as the following optimization problem: In this approach, âj is assigned to a i if C ij has the maximum value in the j th column of C. This straightforward approach has two drawbacks. On one hand, if two or more maximum values occurred in the same row, a reasonable assignment could not be concluded. This happens for instance in far-field antenna array processing when sources are angularly close, in the presence of noise [START_REF] Sahnoun | Joint source estimation and localization[END_REF]. On the other hand, the delivered index is not reliable, since, even if the index is zero, one cannot guarantee A = AP σ Λ. The following toy numerical example illustrates this problem.

argmin σ 1 2 N n=1 ∥a n -âσ(n) ∥ 2 2 = argmax σ N n=1 |a H n âσ(n) |. (4.2) Let C ij = |a H i âj |,
Assume that in an experiment matrix C is:

C =     0.8 0.3 0.1 0.85 0.9 0.5 0.5 0.2 0.7     . (4.3)
The concluded assignment by this method is (â 1 , a 2 ), (â 2 , a 2 ), (â 3 , a 3 ), which is obviously not acceptable because column a 2 is selected twice. Computing the square error via

1 2 3 n=1 ∥a n -âσ(n) ∥ 2
2 by considering the assumption of normalized a n and âσ(n) with respect to L 2 norm and by substituting the values of |a H i âj | from C ij , one obtains 3 -0.85 -0.9 -0.7 = 0.55, which is less than the exact error, 3 -0.8 -0.9 -0.7 = 0.60 (the exact error is given in Section 4.3.3 with the optimal permutation). This example shows that this algorithm outputs a matrix P p that may not be a permutation. This index is always optimistic since it searches in a set of assignments larger than Perm(N ). In fact, if a set A contains a set B, i.e. B ⊆ A, then

min x∈A f (x) ≤ min x∈B f (x) (4.4)
for any function f (x). Therefore, the reported error is always smaller than or equal to the exact error based on the optimal assignment.

In order to avoid a non-acceptable assignment, after detecting the maximum value of each column of C, its row and column can be removed for the rest of the algorithm. In other words, if in j th column of matrix C, C ij is the maximum value, then the i th row and j th column of C will be ignored in the search of the next maximum value.

This is a greedy approach, since the index depends on the order of choosing the maximum values. For example, if this greedy algorithm is applied on matrix C expressed in (4.3), the resulted assignment will be (â 1 , a 2 ), (â 2 , a 1 ), (â 3 , a 3 ) provided that the columns are swept from left to right. However, if the columns are swept in the opposite way, the assignment will be (â 1 , a 1 ), (â 2 , a 2 ), (â 3 , a 3 ). Compared to the optimistic index, the error output by this greedy approach by sweeping from left to right, 3 -0.85 -0.3 -0.7 = 1.15, is larger than the exact error, 3 -0.8 -0.9 -0.7 = 0.60, while by sweeping from right to left, the reported error equals to the exact error 0.6.

By imposing a column ordering, this greedy approach searches a set of assignments smaller than Perm(N ): following (4.4), one can conclude that the error measurement is always pessimistic. Therefore, the reported error is always larger than or equal to the exact error based on the optimal assignment.

Score measure [BBK18]

This index, which is also known as congruence [START_REF] Stegeman | Using the simultaneous generalized schur decomposition as a candecomp/parafac algorithm for illconditioned data[END_REF], is customized for tensors and is applied to evaluate the performance of a tensor decomposition in terms of estimating all the loading matrices (defined below) together. Let us explain the permutation ambiguity by means of a tensor decomposition example called Canonical Polyadic (CP) [START_REF] Comon | Tensors: a brief introduction[END_REF]. The CP decomposition of a third order tensor of rank 2 admits the following form:

T I×J×K = 2 r=1 a (1) r ⊗ a (2) r ⊗ a (3) r , (4.5)
where ⊗ denotes the outer (tensor) product, and a

(1)

r , a (2)
r and a

(3) r are some vectors of size I, J and K, respectively. Equation (4.5) can be represented in a compact form as T = A (1) , A (2) , A (3) , where

A (i) = [a (i) 1 , a (i)
2 ] is called the mode-i loading matrix of T . Observe that the permuted version of loading matrices, i.e.

A (i) p = [a (i) 2 , a (i)
1 ], i = 1, 2, 3, results in the same tensor as T in (4.5).

The score measure of the tensor

T = A (1) , A (2) , A (3) is calculated based on the correlation matrix C = C (1) C (2) C (3)
, where is the Hadamard product (element-wise product) and

C (k) ij ≜ |a (k) i H â(k) j |, k = 1, 2, 3
. This index is also greedy, since the assignment is concluded based on the maximum values of C, which have been chosen in a way explained in Section 4.3.1.2, and the corresponding score is an average of these selected values.

Methods based on graph matching

The optimal assignment (or optimal transport) problem is an old, wellknown and fundamental combinatorial optimization problem [Gal83,Mun57, PC + 19]. The first polynomial time algorithm for optimal assignment problems is the "Hungarian method" [Kuh55] also known as "Kuhn-Munkres" [Mun57, TK04], and the complexity of the algorithm is approximately N 4 flops [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]. This algorithm has been employed in [START_REF] Tichavsky | Optimal pairing of signal components separated by blind techniques[END_REF] for an optimal pairing of the sources in BSS.

The optimal assignment problem can also be considered as a special case of Maximum Weighted Matching (MWM), which is a well-known problem in graph theory, for which several polynomial time algorithms exist [START_REF] Galil | Efficient algorithms for finding maximal matching in graphs[END_REF].

The best exact [DP14, DK69, D + 59] and approximate [START_REF] Harold | Faster scaling algorithms for network problems[END_REF] MWM algorithms cost approximately 8N 3 and N 2 flops, respectively.

Methods based on optimal permutation

Searching for the optimal permutation σ, i.e. for the optimal permutation matrix P ⋆ , described at the beginning of Section 4.3, can be viewed as finding some entries of C such that no pair among them lies in the same row or column, while the sum of these entries is maximum. One can formulate this as the following optimization problem [PC + 19]:

P ⋆ = argmin P p∈R N ×N + i,j D ij P p ij s.t. P p 1 N = P T p 1 N = 1 N , (4.6)
where

D = -C, 1 N is
p ⋆ = argmin p∈R N 2 + d H p s.t. Qp = 1 2N , (4.7)
where

Q = [1 T N ⊠ I N , I N ⊠ 1 T N ] T ∈ R 2N
×N 2 , I N and ⊠ denote the identity matrix of size N and the Kronecker product, respectively. Yet, from Birkhoff's Theorem, the set of bistochastic matrices is a polyhedron1 whose vertices are permutations [HJ99, Theorem 8.7.1]. On the other hand, a fundamental theorem of linear programming [BT97, Theorem 2.7] states that the minimum of a linear objective in a non-empty polytope (i.e. a finite polyhedron) is reached at a vertex of the polytope. This permits to relax the search for a permutation into (4.6) or (4.7): in fact, looking for the best bistochastic matrix will eventually yield a permutation and this is the whole power of this method, which employs linear programming to estimate the permutation.

For example, by employing MWM or the linear program described above, the optimal permutation in experiment (4.3) is the identity matrix.

Comon index [Com94]

Comon's index is a combination of L 1 and L 2 norms, and is calculated as:

ϵ 1 (S) = N i=1 N j=1 |S ij | -1 2 + N j=1 N i=1 |S ij | -1 2 , N i=1 N j=1 |S ij | 2 -1 + N j=1 N i=1 |S ij | 2 -1 .
In [START_REF] Comon | Independent component analysis, a new concept? Signal Proc[END_REF], it has been proved that ϵ 1 is invariant to permutation, i.e. ϵ 1 (A, A) = ϵ 1 (A, AP p Λ). Moreover, it has been shown that ϵ 1 (A, A) = 0 if and only if A = AP σ Λ, where σ is the optimal permutation. However, ϵ 1 can increase enormously, depending on the values of matrix S, hence, this index is not bounded from above.

Moreau-Macchi index [MM94]

The index proposed in [START_REF] Moreau | A one stage self-adaptive algorithm for source separation[END_REF] measures a gap between matrix S and a permutation matrix. It is defined as:

ϵ 2 (S) = N i=1   N j=1 |S ij | 2 (max k |S ik |) 2 -1   + N j=1 N i=1 |S ij | 2 (max k |S kj |) 2 -1 .
Dividing by the maximum value (e.g. (max k |S ik |) 2 ) provides an upper bound for ϵ 2 unlike ϵ 1 . + 96] This performance index takes the form:

Amari index [ACY

ϵ 3 (S) = N i=1   N j=1 |S ij | max k |S ik | -1   + N j=1 N i=1 |S ij | max k |S kj | -1 .
The only difference between Amari and Moreau-Macchi index is the power 2, which exists in ϵ 2 . Therefore, calculating ϵ 3 is less costly compared to ϵ 2 .

In addition, as for ϵ 2 , the division by the maximum value (e.g.

max k |S ik |)
provides an upper bound for ϵ 3 .

CHAPTER 4. PERFORMANCE INDEX

An accurate investigation of indices reviewed in this section reveals that ϵ 1 is not bounded from above. Furthermore, the upper bounds on ϵ 2 and ϵ 3 have not been studied in [MM94, ACY + 96], so that their upper bound cannot be easily interpreted. Even if one normalizes the Amari and Moreau-Macchi indices, the resulted upper bounds are reached when A has equivalent columns and A is identity matrix. Therefore, the upper bounds of the Amari and Moreau-Macchi indices do not correspond to the largest possible angular gap between A and A.

In order to obtain interpretable upper bounds and to reduce computational cost, one may think of replacing S = A † A by C, but in this case the property that ϵ i = 0 is equivalent to ϵ 0 = 0, for i ∈ {1, 2, 3} does not hold anymore.

Our proposed index: CorrIndex

In this section, we introduce, "CorrIndex", which is based on a correlation matrix. Reminding that we define C = |A H A|, where A ∈ C M ×N and A ∈ C M ×N and modulus is applied entrywise. In addition, we assume that the columns of A and A are normalized by their L 2 norms.

Basically, if ϵ 0 (A, A) = 0, N entries of C are one, since |a n | = |â n | and the columns of A and A are normalized to unit L 2 norms. Remember that it is desired that a performance index is zero if and only if ϵ 0 (A, A) = 0.

In order to satisfy these basic requirements in the matrix case, i.e. M > 1,

CorrIndex is defined as follows:

CorrIndex(C) = 1 2N N i=1 | max k C ik -1| + N j=1 | max k C kj -1| . (4.8)
The Remark: It can be also observed that CorrIndex is bounded:

0 ≤ CorrIndex ≤ 1.
According to (4.8), unlike ϵ 2 and ϵ 3 , the upper bound of CorrIndex is easier to interpret when M > 1, since it is achieved when entries of C are minimal (i.e. the largest possible angular distance between A and A). In particular, when M ≥ 2N , C = 0 when all the columns of A and A are orthogonal to each other, which yields CorrIndex = 1. Next, as proved below, the zero lower bound is meaningful, since it corresponds to

ϵ 0 = 0.
The one-row case: On the other hand in the row vector case, i.e. M = 1, as CorrIndex is based on (4.2), we should return back to the basic minimization of finding optimal assignment (σ), which is a restatement of (4.2) as follows: Comparing (4.8) and (4.10) reveals that "max" and "1" have been replaced with "min" and "0", respectively, which helps benefit from the same properties as (4.8) in the vector case.

argmin σ 1 2 N n=1 (a n -âσ(n) ) 2 . ( 4
In the following, it is shown that CorrIndex is invariant to scaling and

permutation, i.e. CorrIndex(A, A) = CorrIndex(A, AP p Λ). Moreover, it is shown that CorrIndex(C) = 0 if and only if A = AP σ Λ, i.e. ϵ 0 (A, A) = 0.
Proposition 1. CorrIndex is invariant to permutation and scaling: As mentioned before, it is hard to assess the relative error made on loading matrices in tensor decompositions, because of scaling and permutation ambiguities [START_REF] Comon | Tensors: a brief introduction[END_REF]. So as to overcome these ambiguities, we can use Cor-rIndex as a performance on estimating loading matrices. However, if we report CorrIndex on each loading matrix separately, it would be an optimistic index, since implicitly a different permutation would be permitted for each loading matrix. Note that according to (2.6), the corresponding columns of loading matrices construct a rank-1 tensor, hence to keep these rank-1 components unchanged, it is required to permute all loading matrices by the same permutation matrix. By permitting a different permutation for each loading matrix, the overall tensor reconstruction error could be less than the case in which all loading matrices have the same permutation.

CorrIndex(A, A) = CorrIndex(AP p Λ, A) = CorrIndex(A, AP p Λ). (4.11) Proof. Assume that C 1 = |A H A| and C 2 = |(AP p Λ) H A|.
CorrIndex(A, A) = 0 ⇐⇒ A = AP σ Λ. (4.12) Proof. Firstly, if A = AP σ Λ, then max k C ik = 1,
In order to have a more reliable performance index, we can apply Cor-rIndex to a matrix X, built upon loading matrices stacked one below the other. In other words, for the tensor described in (4.5), CorrIndex(X, X) applies to:

X =     A (1) A (2) A (3)     , X =      A (1) A (2) A (3)      .

Discussion and computer results

A multi-aspect comparison between CorrIndex and other reviewed methods has been carried out, and is reported in Table 4.1, where the methods of Section 4.3.1 and 4.3.2 are referred by "Greedy" and "Graph", respectively.

The number of multiplications of each stage, i.e. computing the input matrix (C = |A H A| or S = A † A), estimating the permutation and computing the index, are reported. In addition, the last column of Table 4.1 ("Significance of upper bound") indicates if the upper bound makes sense, i.e. returning the maximum index value for the largest distance between A and A. According to Table 4.1, it is inferred that CorrIndex has the lowest computational complexity compared to the others in terms of the number of multiplications besides its theoretical guarantee, its invariance to permutation and scaling ambiguity and its meaningful bounds.

In the rest of this section, we report either the relative error (for the greedy and graph-based methods), which is defined in (4.13) or the indices ϵ i (Moreau-Macchi, Amari, . . . ). As greedy and graph-based methods estimate the permutation, we report the relative error between normalized A and A, by means of estimated matrix P p , as follows:

relative error = ∥A -AP p ∥ F ∥A∥ F , (4.13)
where ∥.∥ F denotes the Frobenius norm. Reminding that P p is aimed to be a permutation matrix, but that it might not be, and the inferred assignment of the columns by P p does not make sense (see Section 4.3.1.1).

The index and computation time of each index in a numerical experiment is reported in Table 4.2 to evaluate the methods practically. This experiment

Greedy N 2 M 0 2M N No Graph N 2 M 8N 3 2M N Yes Linprog N 2 M N 4 2M N Yes Comon [Com94] 11N 3 - 2N 2 No Moreau-Macchi [MM94] 11N 3 - 2N 2 No Amari [ACY + 96] 11N 3 - 2N No CorrIndex N 2 M - 1 Yes
is executed on a laptop with a processor of 3.1 GHz Intel Core i5, 16 GB RAM, running macOS Mojave and MATLAB 2019a.

In order to show the drawbacks of greedy methods, this experiment is done on some matrices, A ∈ R M ×N , whose columns are highly correlated.

For this purpose, a correlation matrix, R N ×N , of the columns of A is designed such that its diagonal and off-diagonal entries are 1 and γ, respectively, where γ is an arbitrary mutual coherence constant among the columns of A. Then, by considering the Cholesky decomposition of R, i.e. R = L T L, and a random orthogonal matrix U M ×N (U can be obtained by the QR decomposition of a random matrix), we set A = U L. In this way, we have a matrix A with the columns having the correlation of γ.

A is generated by permuting randomly the columns of A and adding a noise matrix, W , of the same size as A with i.i.d. entries of Gaussian distribution with zero mean and unit variance, and weighted by the parameter δ. The variance δ 2 of the additive noise is adjusted such that we reach a desired Signal to Noise Ratio (SNR) defined as:

SNR = 10 log 10 i,j A(i, j) 2 i,j δ 2 W (i, j) 2 .
(4.14)

At the end, the columns of A and A are normalized. In the experiment of Table 4.2, M = 150, N = 100, with the mutual coherence constant γ = 0.75, δ = 0.1 (which is equivalent to SNR = -1.76 dB), and U is an orthogonal matrix obtained by concatenating the first N leftsingular vectors of a random matrix whose entries are chosen randomly from a uniform distribution on (0, 1). The reported values are averaged over 50 realizations.

The indices obtained by greedy methods and reported in Table 4.2 explicitly show the effect of coherence of input matrix on these types of methods.

For instance, according to the performed experiment, greedy methods either report less (37%) or larger (105%) error than the exact index (86%).

Note that as greedy methods try first to estimate permutation P p , and then calculate the error between matrices AP p and A; hence, the indices computed by greedy methods may be compared to the exact error computed by graph-based methods. However, comparing other indices such as CorrIndex with the exact error does not make sense, since these indices are intrinsically As it can be interpreted from the reported indices in Table 4.3, the greedy method of [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF][START_REF] Coloigner | Line search and trust region strategies for canonical decomposition of semi-nonnegative semi-symmetric 3rd order tensors[END_REF] does not report zero relative error between A and A, which is not correct. However, in spite of highly correlated columns of A, all indices (except the greedy method of [CLDA09, CKAC14]) demonstrate zero distance between A and A, which is true. Comparing Table 4.2 and Table 4.3 reveals that the greedy methods are much more sensitive to the correlation of the columns than other indices.

CorrIndex is based on (4.2), which tries to minimize the least square error between A and A. Therefore, if the distance between A and A increases due to the additive noise in A, CorrIndex will return a larger value. To show this fact in practice, we performed an experiment whose result is depicted in Fig. 4.1. Generating a random matrix A of dimension 6 × 4, A is obtained by permuting its columns and by adding a noise matrix, W , of the same size as A with independent and identically distributed (i.i.d.) entries of Gaussian distribution with zero mean and unit variance, and weighted by the parameter δ. The variance δ 2 of the additive noise is adjusted such that we reach a desired SNR as described in (4.14). ) is optimistic (resp. pessimistic), since its reported error is smaller (resp. larger) than the exact error. In addition, as SNR increases, this error gets larger, which demonstrates that by decreasing the additive noise, the influence of mutual coherence becomes more effective on the result of greedy methods.

Conclusion

In this chapter, the problem of computing the distance between two matrices up to permutation and scaling ambiguities is addressed. This problem occurs

Introduction

In this chapter, we study constrained and non-constrained tensor decompositions, as well as introducing our own constrained tensor decomposition, Simple Forward-Backward Splitting (SFBS). As mentioned in Chapter 3, the main step in data mining by tensors is to decompose properly the tensor of third order moments. As we will show by simulation in this chapter, constrained tensor decomposition results in a more accurate and more reasonable (in the sense of being in the simplex set) estimation of probabilities of hidden and multi-view variables, which is the first step for some targeted data mining tasks such as unsupervised clustering. In addition, we review some non-constrained decompositions, which have been employed for this purpose, and will compare them experimentally and theoretically with constrained tensor decompositions.

Depending on the application, it is preferred to add some constraints to the tensor decomposition, which results normally in much more accurate and reasonable solutions. Non-negativity, belonging to the simplex set, orthogonality and sparsity are some examples of key constraints often imposed in some applications such as medical image and signal processing [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF],

probability estimation in topic modeling (cf. Section 3.2) and dictionary learning [START_REF] Hsieh | 2d sparse dictionary learning via tensor decomposition[END_REF].

Generally, the algorithms of constrained tensor decomposition are inspired from constrained matrix decomposition (factorization). For instance, the algorithms mentioned in [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF] for Non-negative Tensor Factorization (NTF) are extensions of Non-negative Matrix Factorization (NMF).

Moreover, many of the existing algorithms are based on Alternating Optimization (AO) [SDLF + 17] or its special case, Alternating Least Squares (ALS) [CMDL + 15], in which the data fidelity term in AO is the least square error.

Over the past decade, some algorithms have been proposed for constrained tensor decomposition based on AO or ALS [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF][START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF] Vu | A proximal approach for nonnegative tensor decomposition[END_REF],

where in each step, a constrained minimization over one parameter is carried out. In each step, Alternating Direction Method of Multipliers (ADMM)

[BPC + 11] or proximal methods [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF] have been applied to solve the constrained minimization.

As our proposed method, SFBS, is based on the proximal concept, we limit ourself to review briefly the state-of-the-art algorithms utilizing proximal approaches. Therefore, in this chapter, after reviewing some preliminaries on the proximal concept in Section 5.2, some state-of-the-art algorithms for constrained tensor decomposition based on the proximal concept are reviewed in Section 5.3. In Section 5.5, we propose our algorithm, i.e. SFBS, and provide its convergence analysis. In Section 5.6, several experiments have been performed to compare our algorithm with constrained and unconstrained algorithms, in terms of performance. We also evaluate tensor decomposition methods on a real text data set, called "20 Newsgroups".

In Section 5.7, we discuss briefly about advantages and drawbacks of each method based on the simulations and experiments in Section 5.6. Finally, Section 5.8 concludes the chapter along with explaining some future perspectives.

Proximal concept and approach

In this section, we explain the proximal concept, operators and methods such as Forward-Backward Splitting, which is used in our proposed method in Section 5.5. Firstly, some required definitions and properties of functions are reviewed as follows: For a function f : R n → R and a point x for which f (x) is finite, a 1) φ(0) = 0,

* Lower semi-continuity (lsc) [RW09, Definition 1.5]: Suppose R = [-∞, +∞]. The function f : R n → R is called lower semi- continuous (lsc) on R n if lim ϵ→0 ϵ>0 [ inf x∈B(x,ϵ) f (x)] = f (x) (5
vector v ∈ R n is said to be the sub-gradient of f at x, i.e. v ∈ ∂f (x), if f (x) ≥ f (x) + ⟨v, x -x⟩ + O(|x -x|),
2) φ is differentiable on (0, η),

3) φ ′ (y) ≥ 0 for all y ∈ (0, η),

4) The Kurdyka-Lojasiewicz inequality,

φ ′ (f (x) -f (x * ))dist(0, ∂f (x)) ≥ 1, holds for all x ∈ U ∩ {x|f (x * ) < f (x) < f (x * ) + η}, where dist(.)
denotes the distance function. In the rest of the report, this property is referred to as "KL", in short.

In the context of optimization, the KL property is important, since many problems include functions satisfying this property. As it is explained in Remark 2, many semi-algebraic (cf. Remark 2 for definition of semi-algebraic) functions have the KL property, however,

) def = {(x, f (x)) : x ∈ dom(f )}.
In the rest of this section, the proximal concept and approaches are reviewed.

Proximity operator

The projection of a vector x ∈ R N onto a closed convex set S ⊂ R N is a classical problem in signal processing, which is depicted in Fig. 5.2 and can be formulated as [CZ + 97, Com93, YW82]:

proj S (x) = argmin y∈R N i S (y) + 1 2 ∥x -y∥ 2 2 , (5.3)
where i S is the indicator function defined by:

i S (y) ≜    0 if y ∈ S ∞ if y ̸ ∈ S (
argmin y∈R N f (y) + 1 2 ∥x -y∥ 2 2 , (5.5)
is defined as the proximity operator of the function f ∈ Γ 0 (R N ), and it is denoted by prox f (x). The term 1 2 ∥x -y∥ 2 2 is also called proximal regularization in the literature [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF][START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF]. Thus, the proximity operator of f is prox f : R N → R N , and it is characterized by:

p = prox f (x) ⇔ (x -p) ∈ ∂f (p), ∀(x, p) ∈ R N × R N .
Note that ∂f (p) is replaced by ∇f (p) for differentiable f . The above definition indicates that prox f (x) is a point that minimizes f and simultaneously is as close as possible to x. Therefore, prox f (x) is also called a proximal point of x with respect to f [START_REF] Parikh | Proximal algorithms[END_REF]. See Table 10.2 in [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF] for a list of popular functions and their corresponding proximity operators.

However for the sake of clarity, Table 5.1 provides proximity operators of some functions, which are discussed in this thesis.

Proximal methods

Although proximal methods have first appeared in the work of Martinet [START_REF] Martinet | Régularisation d'inéquations variationnelles par approximations successives[END_REF] in 1970, they have only been utilized in signal processing since 2001 [START_REF] Patrick | Signal recovery by proximal forward-backward splitting[END_REF], and their applications in various fields are getting more and more prevalent. As mentioned before, a proximity operator can be viewed as a generalized projection [START_REF] Parikh | Proximal algorithms[END_REF]. In addition to projection, there are several interpretations, some of which we review from [PB14, CP11]:

-The approximation of gradient step:

If f is twice differentiable at x and its Hessian matrix is positive definite, then it can be shown that as λ → 0:

prox λf (x) = x -λ∇f (x) + O(λ).
Therefore, for small λ, prox λf (x) converges to a gradient step in f with step length λ. Thus, the proximity operator can be interpreted as an approximation of a gradient step for minimizing f provided that λ is small enough. This fact shows a connection between proximity operators and gradient methods, and consequently the proximal operator could be useful in optimization. It can be also expected that λ will have a similar effect as a step size in a gradient method.

-Fixed-point of prox λf :

It can be proved that the fixed points of the proximity operator of f are exactly the minimizers of f [PB14, sec 2.3], i.e. prox λf (x * ) = x * if and only if x * minimizes f . This also reveals a close connection between proximity operators and fixed point theory. Therefore, fixed points of appropriate operators (proximity operators) can be utilized as a solution for optimization problems. This interpretation will be used in Theorem 1 of Section 5.2.3.

-Denoising:

A proximity operator can also be interpreted as a solution for denoising in signal processing applications. Consider a noisy observation y ∈ R N of a signal x ∈ R N , y = x + w, where w models the additive noise.

The best recovery of the original signal x is usually formulated as:

argmin x∈R N f (x) + 1 2 ∥y -x∥ 2 2 , (5.6)
where 1 2 ∥y -x∥ 2 2 is the data fidelity term and f (x) takes into account a prior knowledge about x. According to the definition of proximity operator, (5.6) is exactly equal to prox f (y).

Forward-Backward Splitting

In many signal processing applications, the cost function to be minimized is the sum of two functions where one of them is usually non-differentiable or even non-convex. By following a proximal approach, these kinds of problems Theorem 1 (Forward-Backward Splitting [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]). Suppose f : R N → R ∪ {+∞} is a proper3 lower semi-continuous function which has the KL property and is bounded from below. If f can be split into two parts as f = h+g, where g is lower semi-continuous and h : R N → R is a finite valued, differentiable function with a β-Lipschitz continuous gradient, i.e., ∃β such that:

∥∇h(x) -∇h(y)∥ 2 ≤ β∥x -y∥ 2 ,
then it can be shown [START_REF] Patrick | Signal recovery by proximal forward-backward splitting[END_REF] that the minimizer of f satisfies the following fixed point equation:

x = prox γg (x -γ∇h(x)), (5.7)
where γ ∈ (0, +∞).

Equation (5.7) suggests an iterative approach, called the Forward-Backward Splitting algorithm:

x k+1 = prox γg (x k -γ k ∇h(x k )), (5.8)
where the values of γ k should be chosen from a suitable bounded interval.

Forward-Backward Splitting is a combination of two basic methods: a proximal algorithm and a gradient approach [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]. Actually, if g = 0, then (5.8) is transformed to the gradient method of the differentiable function h. On the other hand, when h = 0, (5.8) is the proximity operator of g, and the minimizer of f is the fixed point of this operator.

Several variations of implementing Forward-Backward Splitting exist and are reported in [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]. Two of them are restated (Algorithm 2 and Algorithm 3) to which we will refer in the rest of this thesis.

In Algorithm 2, several parameters are required to be set by user, such 2: for k = 0, 1, 2, . . . do 3:

γ k ∈ [ϵ, 2 β -ϵ]
4:

y k = x k -γ k ∇h(x k ) 5: α k ∈ [ϵ, 1]
6:

x k+1 = x k + α k (prox γ k g (y k ) -x k )
7: end for implements the iterative approach described in (5.8), except doing an extrapolation with α k in line 6 between the new point obtained from applying the proximity operator and the previous estimation point (x k ).

Algorithm 3 is based on Fast Iterative Shrinkage Thresholding Algorithm (FISTA) proposed in [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF][START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF] and can be considered as a proximal gradient algorithm. In fact, instead of shrinkage function in FISTA, the proximity operator is used.

The main difference between Algorithm 2 and Algorithm 3 is that in Algorithm 3 the proximity operator is not employed on the previous iterate

x k but rather on a very specific linear combination of the previous two iterates, x k and x k-1 [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF].

Although Algorithm 3 is more user-friendly than Algorithm 2 in terms of the number of required parameters to be set, the computation complexity of Algorithm 2 is less than Algorithm 3 due to not calculating some coefficients like t k and λ k in lines 5 and 6 of Algorithm 3.

In many applications (including SFBS in Section 5.5), the function g is an indicator function of a particular set (S), i S , and its proximity operator is a projection onto that set [CP11, Table 10.2]. If the desired set is nonconvex, the projection onto it may not result to a unique point. It has been proved [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF] that in spite of the multi-valued projection, the convergence

y k = z k -β -1 ∇h(z k ) 4:
x k+1 = prox β -1 g (y k ) 5:

t k+1 = 1 + 4t 2 k + 1 2 6: λ k = 1 + t k -1 t k+1
7:

z k+1 = x k + λ k (x k+1 -x k )
8: end for property of Theorem 1 is not influenced. Note that this interesting conclusion is valid only if the assumptions of Theorem 1 are satisfied, the most important being the KL property satisfied by h + i S .

Constrained algorithms based on the proximal concept

In this section, some algorithms that use proximal approaches for constrained (mostly non-negativity) CP decomposition are reviewed. Our proposed method (SFBS in Section 5.5) will be compared to all of them in Section 5.6.1.

Problem formulation

The problem of N order constrained CP decomposition can be expressed as the following minimization:

min λ,A (1) ,...,A (N ) 1 2 ∥T -λ; A (1) , A (2) , . . . , A (N ) ∥ 2 F s.t. C λ , C A (n) , 1 ≤ n ≤ N
where T is a N order tensor of dimension n 1 ×n 2 ×. . .×n N and rank R with the coefficient vector λ R×1 and N loading matrices A

(1)

n 1 ×R , A (2)
n 2 ×R , . . . , A

n N ×R . In addition, C λ , C A (n) represent the constraint over the coefficient vector λ and loading matrix A (n) , respectively.

Since minimizing the cost function of the CP decomposition over all loading matrices is a non-convex problem, a common strategy to transform it to a sequence of convex problems (if constraints are already convex) is Alternating Optimization (AO) framework [SDLF + 17] or Block Coordinate Descent (BCD) [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]. In AO, by fixing all the loading matrices (by initialization or using their previous estimation) except one of them, one tries to minimize the cost function over just one loading matrix.

For example, ALS steps for a third order constrained tensor decomposition of dimensions n 1 ×n 2 ×n 3 and rank R, T = 1 R ; A (1)

n 1 ×R , A (2) n 2 ×R , A (3)
n 3 ×R , in an unfolding form, are as follows:

A (1) k+1 = argmin A n 1 ×R 1 2 ∥T (1) n 1 ×n 2 n 3 -A(A (3) k ⊙ A (2) k ) T ∥ 2 F s.t. C A (A) A (2) k+1 = argmin B n 2 ×R 1 2 ∥T (2) n 2 ×n 1 n 3 -B(A (3) k ⊙ A (1) k+1 ) T ∥ 2 F s.t. C B (B) A (3) k+1 = argmin C n 3 ×R 1 2 ∥T (3) n 3 ×n 1 n 2 -C(A (2) k+1 ⊙ A (1) k+1 ) T ∥ 2 F s.t. C C (C) (5.9)
where C • (•) represents the desired constraint over loading matrices, and it can be added to its corresponding fidelity term in the form of generalized regularization, e.g. an indicator function as defined in (5.4).

Non-negativity is a relevant constraint that has been employed in the literature, but in this thesis, we discuss also other useful constraints such as simplex set and sparsity.

Alternating Optimization-Alternating Direction Method of Multipliers (AO-ADMM) [HSL16]

As mentioned in [PB14, CP11], Alternating Direction Method of Multipliers (ADMM) [BPC + 11] can be considered as a special case of the proximal method. In [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF], a constrained CP decomposition by means of ADMM is discussed.

AO-ADMM attempts to minimize each step of (5.9) by ADMM. The details of AO-ADMM for constrained CP decomposition can be found in [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF].

The proposed algorithm in [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] is capable of applying several constraints on loading matrices, such as non-negativity, sparsity, smoothness, cardinality, etc.

The convergence of AO is briefly reviewed in [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] based on the works of [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF][START_REF] Meisam Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF]. Although it is mentioned that adding a proximal regularization (cf. (5.5)) provides a unique solution which serves the convergence of BCD, this regularization is omitted in the formulation of AO-ADMM for the sake of convenience. The convergence of each step which is an ADMM algorithm can be shown for convex functions [HSL16,BPC + 11] and also recently for some non-convex functions such as ℓ q , 0 < q < 1 [WYZ19], but not yet for some others such as ℓ 0 . Note that ℓ q norm of the vector x of dimension n for q > 0 is defined as ( n i=1 |x i | q ) 1 q ; however ℓ 0 (x) pseudo-norm is the number of the non-zero elements of x.

Alternating Proximal Gradient (APG) [XY13b]

In [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF], the convergence of a general BCD algorithm is studied, where for updating each block of unknown variables, three kinds of updates are investigated, namely original, proximal and prox-linear. In two of them, proximal regularization is considered. This general BCD algorithm can be applied on a vast variety of constrained optimization, including constrained tensor decomposition.

Based on the comparison done in [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF], prox-linear updating (i) yields better objective values than original and proximal updating, (ii) is easier to compute and (iii) often allows closed form solutions.

APG does not explicitly follow Forward-Backward Splitting procedure, nevertheless its algorithm for non-negativity constraint ends in projection to non-negative orthant along with updating a particular coefficient as in line 5 of Algorithm 3.

The convergence analysis in [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF] reveals that BCD with original updating requires strong convexity to provide a particular convergence guarantee4 , while by proximal and prox-linear updating, the global convergence to a critical point is achieved under mild conditions such as satisfying KL inequality or Lipschitz continuity. As APG, FastNTF-APG utilizes prox-linear updating, and its idea to speed up convergence can be explained by following each step of (5.9). Without loss of generality, let us explain the idea motivating the first step. In order to reduce the computational complexity, matrix T (1) in T (1) (A

Fast

(3) k ⊙ A (2)
k ) T is replaced with its low rank approximation. This is also efficient for filtering the noise out of T (1) .

As in APG, it is also needed in FastNTF-APG to compute the Lipschitz constant of the gradient of 1 2 ∥T (1) -A(A

k ⊙A

(2) k ) T ∥ 2 F w.r.t. A (considering again the update of the first loading matrix, to fix the ideas). It is mentioned in

[ZZZ + 16] that the Lipschitz constant is ∥(A (3) k ⊙ A (2) k ) T (A (3) k ⊙ A (2)
k )∥ F , but as we prove in Appendix A, the Lipschitz constant is in fact the spectral norm 5 of {(A

(3) k ⊙ A (2) k ) T (A (3) k ⊙ A (2) k )}.
Let us remark a point about the importance of the Lipschitz constant of the gradient of the fidelity terms in (5.9). To employ the Forward-Backward Splitting algorithm for minimizing the objective functions in (5.9), the fidelity term (e.g. 1 2 ∥T (1) -A(A

(3)

k ⊙ A (2) k ) T ∥ 2 F )
roles as the function "h" in Theorem 1. In Theorem 1, "β" is the Lipschitz constant of the gradient of "h". On the other hand, "γ" is the coefficient of ∇h in the fixed point equation (5.7) in the Forward-Backward Splitting algorithm (cf. Theorem 1). As we shall see in Section 5.5.4 (cf. Theorem 2), by choosing γ ∈ [0, 1 β ], the convergence of Theorem 1 is guaranteed. Therefore, computing "β" (or the Lipschitz constant of the gradient of the fidelity term) has a crucial role in the convergence of algorithms, which are based on the Forward-Backward Splitting algorithm.

Block Coordinate Variable Metric Forward-Backward (BC-VMFB) [CPR16, VCTMM17, VCTM + 17]

BC-VMFB consists of two main steps: a gradient step related to the data fidelity, which is assumed to be differentiable and has a β-Lipschitz gradient, and a proximal step linked to the regularization term, for which a new proximity operator should be calculated. This proximity operator of the function φ is associated with a symmetric positive definite matrix Z by the following definition [START_REF] Vu | A proximal approach for nonnegative tensor decomposition[END_REF]:

prox Z,φ (v) = argmin u 1 2 ∥u -v∥ 2 Z + φ(v)
where ∥x∥ 2 Z = ⟨x, Zx⟩, and ⟨•, •⟩ is the inner product. In the above definition of the proximity operator, Z is called preconditioning matrix [START_REF] Vu | A proximal approach for nonnegative tensor decomposition[END_REF]. DECOMPOSITIONS The definition of the proximity operator given in (5.2) may be obtained from the above definition, if Z is the identity matrix. It is observed empirically in [START_REF] Chang | Variational phase retrieval with globally convergent preconditioned proximal algorithm[END_REF] that utilizing preconditioning matrix speeds up the convergence of Proximal Alternating Linearized Minimization (PALM) [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF].

Let us bring two points about non-negative tensor decomposition by means of BC-VMFB. Firstly, the convergence of BC-VMFB is studied in [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF],

and it is shown that if the cost function satisfies the KL property, then its convergence to a critical point is guaranteed. Secondly, the algorithm proposed in [START_REF] Vu | A proximal approach for nonnegative tensor decomposition[END_REF] is designed to handle vectors, so each loading matrix has to be vectorized before utilization.

Non-constrained algorithms for symmetric tensor decompositions

As mentioned in Section 3.3.1.1, in order to estimate the probabilities of hidden variables and the conditional probabilities of multi-view variables (which are denoted by φ and A, respectively), it is required to decompose the third order moments tensor (T ). The focus of this section is on the nonconstrained decomposition for this purpose. Remind that the third order moments (T ) and the second order moments (P ) are related to the desired probabilities according to (3.6) and (3.5), which are worth to be mentioned again:

P = K k=1 φ k a k ⊗ a k T = K k=1 φ k a k ⊗ a k ⊗ a k .

Robust tensor power method [AGH + 14]

We describe in this section the approach of [AGH + 14], called "Robust tensor power method", whose drawbacks we shall investigate experimentally in Sec-tion 5.6.2. This method of decomposition is dedicated to symmetric tensors and to the problem described in Section 3.3, where the empirical estimation of the second order moments is also available in addition to the tensor of third order moments. In [AGH + 14], the authors propose to use the two moments (P and T ) defined in (3.3) and (3.4) with exactly the same encoding explained in Section 3.3.1.1. Therefore, these moments satisfy (3.5) and (3.6).

Matrix P is theoretically positive semi-definite, since P is a covariance matrix (note that φ k are non-negative numbers in (3.5)). Hence in a similar way as had been done for Blind Source Separation [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF][START_REF] Comon | Independent component analysis, a new concept? Signal Proc[END_REF], there exists a "whitening" matrix W such that W T P W = I, where I is the identity matrix; W can theoretically be easily obtained from the EigenValue

Decomposition (EVD) P = U DU T , U T U = I, by setting W = U D -1/2 . Therefore, K k=1 ãk ãT k = I, in which ãk def = √ φ k W T a k .
In other words, the matrix A ∈ R K×K containing ãk as its columns is orthogonal, i.e. A A T = I.

As mentioned in Section 3.4, in practice, the matrix P can be estimated via a particular encoding of multi-view variables x ℓ and a moment estimator for approximating (3.3), and consequently, P may have negative eigenvalues (of course it is not the case with the encoding explained in Section 3.3.1.1 and estimators described in Section 3.4) because of estimation errors. This issue was not investigated in [AGH + 14], because they used left singular vectors of P instead of eigenvectors to construct W (cf. Algorithm 4). Nevertheless, the issue still remains, since left and right singular vectors are not imposed to be the same. Next, this whitening matrix is applied to the tensor T to yield:

T def = T • 1 W • 2 W • 3 W ,

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

which have the following element-wise definition according to (2.2) [START_REF] Comon | Tensors: a brief introduction[END_REF]:

T (r, t, s) = r ′ ,t ′ ,s ′ T (r ′ , t ′ , s ′ )W (r, r ′ )W (t, t ′ )W (s, s ′ ).
This new tensor satisfies

T = K k=1 φ -1/2 k ãk ⊗ ãk ⊗ ãk .
The conclusion is that T ideally admits an orthogonal CP decomposition; see e. Unfortunately, except for the tensor power iteration [DLC + 95], the pseudo code in [AGH + 14] is not complete in terms of describing the entire algorithm.

Our description in Algorithm 4 hopefully fills this lack. We used this algorithm in our subsequent computer experiments in Section 5.6.2.

We shall show experimentally in Section 5.6.2 that in spite of its name, Robust tensor power method is not robust to the additive noise. Actually, robustness in [AGH + 14] is considered as robustness to the initial point of eigenvector (cf. line 4 of Algorithm 4). It is proved that any initial vector, close enough to the solution, will converge to an eigenvector of T

Shifted Symmetric Higher-Order Power Method (SS-HOPM) [KM11]

Robust tensor power method [AGH + 14] is based on Symmetric Higher-Order Power Method (S-HOPM), which is a rank-1 approximation and tries to find 6 Recall that there exist several definitions of tensor eigenvectors [START_REF] Lim | Singular values and eigenvalues of tensors: a variational approach[END_REF][START_REF] Qi | Tensor Analysis, Spectral Theory and Special tensors[END_REF]. The definition used in [AGH + 14] -and hence here -is T • v • v = λv, which has the undesirable property that λ depends on the norm of v. In fact, if (λ, v) is an eigenpair, then so is (αλ, αv) for any nonzero α. This is analyzed in e.g. [START_REF] Qi | Tensor Analysis, Spectral Theory and Special tensors[END_REF].

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

Algorithm 4 Robust tensor power method [AGH + 14] Input: The empirical estimation of moments, i.e. T and P

Output:

The estimation of desired probabilities, i.e. A and φ 1:

Whitening: P = U ΣV T ; W = U Σ -1/2 ; T def = T • 1 W • 2 W • 3 W 2: for k = 1, . . . , K do 3:
for m = 1, . . . , M 1 do 4:

1) draw an initial value for vector of unit 2-norm for âk .

5:

2) compute M 2 power iteration updates in norm 2, i.e.:

6: t = T • 2 âk • 3 âk ; φ k ← ∥t∥ 2 ; âk ← t φ k . 7:
end for 8:

Pick the trial (m) having the largest φ k .

9:

Refine âk , φ k by M 2 extraneous iterations.

10:

Deflation: T ← T -φ âk ⊗ âk ⊗ âk 11: end for 12: De-whitening: B = (W T ) † ; a k = φ k B âk ; φ k = 1 φ 2 k for k = 1, . . . , K 13: Back to norm 1: α = ∥a k ∥ 1 ; a k ← a k α ; φ k ← φ k α for k = 1, . . . , K (our suggestion)
the largest eigenvalue of a symmetric tensor [DLC + 95]. Robust tensor power method generalizes S-HOPM to a rank-K approximation by performing deflation after estimating each of K eigenvalues. S-HOPM does not have the guarantee of convergence, unless for an even-order tensor whose corresponding cost function is convex. Shifted Symmetric Higher-Order Power Method (SS-HOPM) improves the convergence of S-HOPM by adding an additive adjustable term, called shift, to the cost function. Although SS-HOPM like S-HOPM requires also the order of tensor to be even for its convergence analysis, it relaxes the convexity assumption on the cost function.

Singular Value based Tensor Decomposition (SVTD) [RCG18]

In spite of other methods described heretofore, Singular Value based Tensor Decomposition (SVTD) is based on matrix decomposition rather than tensor decomposition. In addition, unlike Power method and its variants which utilize the eigenvalues and the eigenvectors of the whitened third order moments, SVTD employs its singular values and singular vectors. Therefore, as Robust tensor power method, SVTD also requires the second order moments for whitening the third order moments. Additionally, the first order moments (or simply the average of multi-view variables) is also needed in SVTD to estimate φ.

SVTD can be viewed as the joint diagonalization of the second order moments (P ) and a particular slice of the whitened third order moments (T ); see [RCG18, Algorithm 1] for much more details. The authors of [START_REF] Ruffini | A new method of moments for latent variable models[END_REF] mentioned that SVTD has some advantages over Robust tensor power method in terms of executing time and required memory space. However, the experiments in [START_REF] Ruffini | A new method of moments for latent variable models[END_REF] reveal that its performance is the same as that of Robust tensor power method.

Proposed Tensor Decomposition Scheme:

Simple Forward-Backward Splitting (SFBS)

In this section, we describe SFBS, the method we propose for constrained (e.g. non-negativity) CP decomposition. First, in Section 5.5.1, the cost function, its solution by means of Forward-Backward Splitting and SFBS algorithm, which is based on Algorithm 2 are introduced. Then, some prevalent constraints such as non-negativity and the simplex space along with required modifications in SFBS algorithm for these particular constraints are discussed in Section 5.5.2. In Section 5.5.4, the convergence theorem for SFBS is stated, whereas some required details of calculations can be found in Appendix A.

Formulation and algorithm

In section 2.3.1, CP decomposition and some constraints are described. Let us now formulate constrained CP decomposition and some prevalent constraints in this section.

Consider the N -th order tensor T ∈ R I 1 ×I 2 ×...I N of rank R. Assume that N ) , where λ ∈ R R + and A (n) ∈ R In×R . A general problem of the constrained CP decomposition of T can be formulated as follows:

T = λ; A (1) , A (2) , . . . , A (
min λ,A (n) 1 2 ∥T -λ; A (1) , A (2) , . . . , A (N ) ∥ 2 F (5.10) s.t. C λ (λ), C A (n) (A (n) ), 1 ≤ n ≤ N,
where C λ (λ), C A (n) (A (n) ) are, respectively the constraints on the vector λ (including the above mentioned constraint, i.e. λ ∈ R R + , such as belonging to the simplex set) and the matrix A (n) .

As mentioned in Section 5.3.2, a common strategy is to solve (5.10) via ALS. Moreover, a constrained optimization can be transformed into an unconstrained one by adding the indicator function of the constraint set to the cost function. To be more precise, at the n th step of ALS for solving (5.10), we have:

min A (n) 1 2 T -A (1) , . . . , A (n) , . . . , A (N ) 2 F + i C A (n) (A (n) ), (5.11) where i C A (n) (A (n)
) is defined as follows:

i C A (n) (A (n) ) =    0 if A (n) ∈ C A (n) ∞ if A (n) ̸ ∈ C A (n)
.

The vector λ is omitted in (5.11), since it can be calculated by normalizing loading matrices (A (n) ). Otherwise, the vector λ as one of the unknown variables can be optimized in one of the steps of ALS.

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

Define W ≜ (A (N ) ⊙ . . . ⊙ A (n+1) ⊙ A (n-1) ⊙ . . . ⊙ A (1) ) T . Then by the mode-n unfolding of (5.11), we have:

min A (n) 1 2 ∥T (n) -A (n) W ∥ 2 F + i C A (n) (A (n) ).
(5.12)

Observe that i C A (n) (A (n)
) is a lower semi-continuous function (cf. the definition of lsc in page 70) for many prevalent constraints such as nonnegativity and the simplex set (see Section 5.5.2 for the definition of these constraints) and

1 2 ∥T (n) -A (n) W ∥ 2
F is finite valued7 , differentiable and β-Lipschitz continuous gradient where β = ∥W W T ∥ σ denotes the spectral norm of matrix W W T (see Appendix A for calculations). As mentioned before, the spectral norm of a matrix is equal to its maximum singular value.

Since W W T is a symmetric matrix, its singular values are the squared of those of W . Moreover, the cost function in (5.12) is proper, lower semicontinuous with the KL property [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF] (see Section 5.5.4 for more explanation). Consequently, all the required assumptions of Theorem 1 are satisfied for (5.12), and according to this theorem, the minimizer of (5.12) is the convergence point of the following fixed point equation:

A (n) = prox γi C A (n) {A (n) -γ(A (n) W W T -W T (n) T )}.
Since the proximity operator of γi C A (n) is the projection onto C A (n) , we have:

A (n) = proj C A (n) {A (n) -γ(A (n) W W T -W T (n) T )}.
(5.13) SFBS is described in Algorithm 5. Although Algorithm 5 is based upon Algorithm 2, we experimentally find the proper values of some parameters of Algorithm 2 such as γ k and α k to obtain the best results. As it is expressed in Algorithm 5, we ignore ϵ and fix the value of γ k = e β in all iterations over k (we experimentally observed that e = 1.4 or e = 1.9 are almost always proper values). In addition, we remove the effect of α g by setting it to 1 in the linear updating of the estimated variable (Line 8 in Algorithm 5).

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

Algorithm 5 The algorithm of SFBS

Input: T , C A (n) , initial A (n) 0 , n ∈ [1, . . . , N ], e Output: Estimated A (n) , n ∈ [1, . . . , N ] 1: repeat 2:
for n = 1, 2, . . . , N do 3:

W = (A (N ) ⊙ . . . ⊙ A (n+1) ⊙ A (n-1) ⊙ . . . ⊙ A (1) ) T 4: β = {max(singular value(W ))} 2 5:
set γ = e β and choose α g .

6:

for g = 0, 1, 2, . . . do 7:

Y = A (n) g -γ(A (n) g (W W T ) -T (n) W T ) 8: A (n) g+1 = A (n) g + α g (proj C A (n) (Y ) -A (n) g ) 9:
end for 10:

end for 11: until some termination criterion

As SFBS, compared to the previous algorithms, requires less parameters to be set (especially compared to APG, which needs a calculation of a coefficient in each iteration), and also it is easy enough to understand and implement, we call it Simple Forward-Backward Splitting (SFBS). Although detailed descriptions about the advantages of SFBS compared to other methods (discussed in this thesis) are provided in Section 5.7, we list them below in brief to complete the current section:

• Unlike other constrained methods, it is explained that how SFBS can handle a variety of constraints such as the simplex set on all loading factors along with a coefficient vector, λ,

• Unlike Robust tensor power method, SFBS is capable of handling overcomplete cases (when the number of hidden variables, K, is larger than the number of multi-view ones, D) in estimating probabilities (cf.

Section 3.3.1),

• Compared to APG, as one of the most efficient methods based on

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

proximal concept, SFBS performs better in noisy scenario in terms of relative reconstruction error as well as estimating loading factors,

• Unlike AO-ADMM, a complete convergence analysis is provided in Section 5.5.4.

Some constraints

As mentioned before, the SFBS is not limited to any particular constraint.

Any constraint whose indicator function satisfies the assumptions of Theorem 1 can be considered. In this section, the required modifications of Algorithm 5 are described, for some widespread constraints, all of which are lsc (cf. the definition of lsc in page 70).

Non-negativity

Non-negativity is one of the most common constraints in the literature.

In many applications such as image processing [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF] or chemometrics [START_REF] Lim | Nonnegative approximations of nonnegative tensors[END_REF], the non-negativity of loading matrices is essential and helpful in the performance of tensor decomposition. The non-negativity constraints of (5.10) are expressed as:

A (n) ∈ R In×R + and λ ∈ R R + .
The projection to the non-negative orthant is done with the max operator, thereby the line 8 of Algorithm 5 would be A

(n) g+1 = A (n) g +α g (max(Y , 0)- A (n) g ).
In other words, the proximity operator of the non-negativity constraint retains the non-negative elements of the array and replaces its negative values with zero.

As mentioned in Section 3.3.1.1, in order to estimate the probabilities of hidden variables and the conditional probabilities of multi-view variables (which are denoted by φ and A, respectively), it is required to decompose the third order moments tensor (T ). One possibility is to apply constrained tensor decompositions such as SFBS. Although it is better to consider φ and A belonging to a simplex set, decomposing the third order moments tensor under the non-negativity constraint of φ and A would be acceptable enough. DECOMPOSITIONS Hence, one can use SFBS with the non-negativity constraint to estimate the desired probabilities in the problem described in Section 3.3.

Simplex set

In the applications involving the probability estimation or the distribution approximation, the simplex set (or the probability simplex) constraint unavoidably appears (cf. Section 3.3.1.1). A vector x ∈ R P belongs to the simplex set S P , if {x T 1 = 1, x i ≥ 0, i ∈ [1, . . . , P ]}. This constraint can be written as {A (n)

(:, j) ∈ S In , j ∈ [1, . . . , R], n ∈ [1, . . . , N ]} and λ ∈ S R .
An algorithm for projecting a vector onto the simplex set is proposed in [START_REF] Condat | Fast projection onto the simplex and the ℓ 1 ball[END_REF]. Therefore, the line 8 of Algorithm 5 would include a projection algorithm to the simplex, whose input is Y .

An implementation trick. As simplex set is a crucial constraint in text mining application (estimating probabilities), here we propose a practical trick to handle this constraint over all parameters (the coefficient vector and all loading matrices).

An efficient way to apply the simplex constraint to all the columns of loading matrices A (n) , n ∈ [1, . . . , N ] and to the coefficient vector λ is to first combine λ with one of the loading matrices, let us say A (N ) . By combination, we mean A (N ) λ = A (N ) Diag(λ), where Diag(λ) is a diagonal matrix containing λ on its diagonal. Then, the simplex constraint may be applied to every column of every matrix A (n) , n ∈ [1, . . . , N -1], and only to the vectorization of matrix A (N ) λ . Each entry λ r of λ is eventually obtained by normalizing the rth column of matrix A (N ) λ with respect to ℓ 1 norm, and the resulting normalized matrix yields an estimation of A (N ) . It can then be proved easily that the estimated λ and every column of loading matrix A (N ) indeed lie in the simplex set.

As mentioned before, constrained tensor decompositions such as SFBS is a reasonable choice to estimate the desired probabilities in the problem described in Section 3.3. Since we seek the probabilities, it is expected that φ and A belong to a simplex set. Following the instructions of ap-

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

plying the simplex set constraint described above, one can utilize SFBS under the simplex set constraint to decompose the third order moments tensor. If we denote the obtained coefficient vector and the loading factors by λ, A (1) , A (2) and A (3) , the estimated φ and A can be considered simply as λ and A (1) + A (2) + A (3) 3 , respectively.

Sparsity with ℓ 0 pseudo-norm (cardinality)

In some applications such as two-dimensional dictionary learning [START_REF] Hsieh | 2d sparse dictionary learning via tensor decomposition[END_REF] and co-clustering [START_REF] Evangelos | Coclustering as multilinear decomposition with sparse latent factors[END_REF], a constraint on the number of the non-zero elements of loading matrices (or cardinality) is needed. This constraint is measured with ℓ 0 pseudo-norm and known as sparsity in compressive sensing [START_REF] David | Compressed sensing[END_REF].

In the penalized form, an ℓ 0 pseudo-norm term is added to the cost function [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Although ℓ 0 is non-convex, the resulted cost function can be solved by means of Theorem 1. As mentioned before, contrary to AO-ADMM, the convergence of algorithms based on Theorem 1, including Algorithm 5, is guaranteed for the ℓ 0 pseudo-norm constraint.

In order to use Algorithm 5 with a cardinality constraint, the proximity operator of ℓ 0 should be calculated. It is mentioned in [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF][START_REF] Parikh | Proximal algorithms[END_REF] that the proximity operator of ℓ 0 is a function called hard-thresholding, defined as follows:

f h η (y) ≜    y |y| > η 0 |y| ≤ η . (5.14)
To be more precise, according to (5.5) the proximity operator of ℓ 0 is defined as: In order to minimize LASSO with Algorithm 5, the proximity operator of ℓ 1 should be computed. In [START_REF] Parikh | Proximal algorithms[END_REF], it is stated that according to (5.5), the

prox γ∥•∥ 0 (x) = argmin y∈R N 1 2 ∥x -y∥ 2 2 + γ∥y∥ 0 , (5.15) which is equal to [f h √ 2γ (x i )] N i=1 (a vector of size N whose elements are f h √ 2γ (x i ) for i ∈ [1, . . . , N ]).
proximity operator of ℓ 1 (x) is f s √ 2γ
, which should be applied element-wise on x. f s η is called soft-thresholding, and is defined as:

f s η (y) ≜        y -η y > η 0 |y| ≤ η y + η y < -η . ( 5

.16)

As for ℓ 0 , f s √ 2γ should be applied on each element of Y in the line 8 of Algorithm 5, instead of a projection.

Constrained CP decomposition based on Proximal Forward-Backward Splitting without AO [NMSC21]

Although all the algorithms discussed so far, including SFBS, are based on AO, some other constrained decompositions such as our work in [START_REF] Nazih | Computation of low-rank tensor approximation under existence constraint via a forward-backward algorithm[END_REF] update all the loading factors together. To be more precise, in AO framework, one tries to update one of the loading factors in each step, while the other loading factors are constant (either by initialization or by updating from previous iterations), the proposed algorithm in [START_REF] Nazih | Computation of low-rank tensor approximation under existence constraint via a forward-backward algorithm[END_REF] tries to estimate all the loading factors, simultaneously. To do this, firstly, all loading factors are concatenated into a one variable matrix, like X defined in (5.22).

Then, Forward-Backward Splitting is applied to minimize a constrained objective over the vectorization of X. For instance, in our work [START_REF] Nazih | Computation of low-rank tensor approximation under existence constraint via a forward-backward algorithm[END_REF], the considered constraint for CP decomposition is the coherence of the columns of loading factors.

Convergence guarantee

In this section, we analyse the convergence of the SFBS algorithm. SFBS, like other methods reviewed in Section 5.3, utilizes BCD or AO (to be exact, ALS). Therefore, the convergence analysis is firstly studied for AO in Section 5.5.4.1, and then, the convergence guarantee of each step (or block) of AO (or BCD) is investigated in Section 5.5.4.2.

The convergence of AO with Proximal Minimization

The overall convergence of AO when the minimization in each step is carried out by means of the proximal concept is proved in [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. Since the convergence analysis provided in [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] is suitable for the algorithm based on Proximal Forward-Backward, we refer the overall convergence of SFBS (which is also based on Proximal Forward-Backward) to the analysis performed in [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], namely, Proximal Alternating Linearized Minimization (PALM).

Let us rewrite (5.10) in an unconstrained manner by means of constraints regularization (i.e. indicator functions or ℓ p norms), and mention some new notations, which help us to discuss the required assumptions of PALM more efficiently:

Ψ(λ; A (1) , A (2) , . . . , A (N ) ) = f d (λ; A (1) , A (2) , . . . , A (N ) )+ N n=1 r n (A (n) )+r λ (λ),
(5.17) where f d (.) is the fidelity term as follows:

f d (λ; A (1) , A (2) , . . . , A (N ) ) = 1 2 ∥T -λ; A (1) , A (2) , . . . , A (N ) ∥ 2 F ,

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

and r n (.) or r λ (λ) could be indicator functions or ℓ p norm, depending on the constraints on A (n) and λ.

In the sequel, the required assumptions of PALM are reviewed, and we shall explain how SFBS, which follows PALM procedure, satisfies these assumptions by considering (5.17). Let us first bring up some required points through our explanation over assumptions.

Based on the definition of KL on page 71, some classes of functions that satisfy the KL property are summarized in [XY13b, Section 2.2], among which we are interested in a particular class of functions, namely semialgebraic. The exact definition of semi-algebraic functions and some of their properties can be found in Remark 2, which will be practical in the explanation over the assumptions of PALM.

Remark 2. A set D ⊂ R d is called semi-algebraic if it can be represented as follows [XY13b, Section 2.2] D = s i=1 t j=1 {x ∈ R d |p ij (x) = 0, q ij (x) > 0}, (5.18)
where p ij , q ij are polynomial functions. A function is semi-algebraic if its graph (see page 72 for definition) is a semi-algebraic set.

Some of the elementary properties of semi-algebraic functions are noted below:

1-Semi-algebraic functions satisfy the KL property.

2-The sum of semi-algebraic functions is semi-algebraic.

3-The sum of a polynomial and a semi-algebraic function is semi-algebraic.

4-The indicator functions of a polyhedral set (such as a non-negative set and the probability simplex) is semi-algebraic.

Assumption 2 (i): inf Ψ(.) > -∞, inf r n (.) > -∞, inf r λ (λ) > -∞: these assump- tions are satisfied, since f d ≥ 0, r n (.) ≥ 0, r λ (λ) ≥ 0.
(ii): f d is globally Lipschitz (see below for a detailed definition): since

f d is quadratic, it is globally Lipschitz.
Globally Lipschitz [BST14, Assumption 2 (ii)]. For any fixed y, the function x → G(x, y) is C 1,1 L 1 (y) , namely the partial gradient ∇ x G(x, y) is globally Lipschitz with moduli L 1 (y), that is

∥∇ x G(x 1 , y) -∇ x G(x 2 , y)∥ ≤ L 1 (y) ∥x 1 -x 2 ∥, ∀x 1 , x 2 ∈ R n Likewise, for any fixed x, the function y → G(x, y) is C 1,1 L 2 (x) .
(iii): The following inequalities are trivially fulfilled, since Assumption 2 (ii) is satisfied (cf. [BST14, Remark 3 (iii)]).

Inequalities 3.5 and 3.6 from [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF].

For i = 1, 2, there exists λ - i , λ + i > 0, such that inf{L 1 (y j ) : j ∈ N} ≥ λ - 1 and inf{L 2 (x j ) : j ∈ N} ≥ λ - 2 sup{L 1 (y j ) : j ∈ N} ≥ λ + 1 and sup{L 2 (x j ) : j ∈ N} ≥ λ + 2 ,
where j is the index denoting a PALM iteration.

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

(iv): The following inequality is satisfied for

f d , whenever f d is a C 2 function (cf. [BST14, Remark 3 (iv)]): since f d is quadratic, it is a C 2 function.
Inequality 3.7 from [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. ∇G is Lipschitz continuous on bounded subsets of R n × R m . In other words, for each bounded subset B 1 × B 2 of R n × R m , there exists Z > 0 such that for all

(x i , y i ) ∈ B 1 × B 2 , i = 1, 2: ∥(∇ x G(x 1 , y 1 ) -∇ x G(x 2 , y 2 ), ∇ y G(x 1 , y 1 ) -∇ y G(x 2 , y 2 ))∥ ≤ Z∥(x 1 -x 2 , y 1 -y 2 )∥.
Assumption 3 Ψ(.) satisfies the KL property: according to Remark 2, for most of prevalent constraints such as non-negativity, the simplex set and ℓ p norm (p = 0, 1, 2, ∞), Ψ(.) is semi-algebraic, hence satisfies the KL. In the rest of this thesis, we call this assumption "KL assumption".

The convergence of each step of AO

The convergence of each step of AO in SFBS, which is Forward-Backward Splitting for each block of variable (cf. Theorem 1), is next discussed. It is usually expected to have the global convergence when the cost function to be minimized is convex [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF], which means that the algorithm generates a converging sequence to the solution regardless of the starting point. However, if the objective function is non-convex, the monotonicity of the sequences generated by descent methods will be broken and oscillatory behaviors may appear [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF]. In order to achieve the convergence in such cases, it is necessary to limit ourselves to functions with some particular properties, such as KL [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF].

The convergence of Forward-Backward Splitting in Theorem 1 is proved in [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF], and the provided convergence analysis is not only applicable for continuous and convex cost functions, but is also usable for non-smooth and non-convex functions (or a non-convex set of constraints). It has been proved

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

that the sequence x k generated by Theorem 1 converges to a critical point of f = h + g, if the mentioned sequence is bounded [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF]. Therefore, the convergence guarantee of each step of SFBS results from the following theorem, quoted from [ABS13]:

Theorem 2 (The convergence of Forward-Backward Splitting [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF]).

Suppose f : R N → R ∪ {+∞} is a proper lower semi-continuous function, which has the KL property and is bounded from below. Assume that f can be spilt into two parts as f = h + g, where g is lower semi-continuous and h : R N → R is a finite valued, differentiable function with a β-Lipschitz continuous gradient.

If the sequence generated by Algorithm 5 is bounded, then this sequence will converge to a critical point of f . In addition, by choosing γ ∈ [0, 1 β ] (γ is the coefficient of ∇h in the fixed point equation (5.7), and β is the Lipschitz constant of ∇h), the values of the cost function are not increasing.

Let us explain that the assumptions in Theorem 2 are completely satisfied for each step of SFBS. According to Theorem 2, the objective function in each step, i.e. (5.12), should be a proper lower semi-continuous (lsc) function, which has the KL property and is bounded from below. As mentioned before, it is proper as its domain is not null. This function is lsc, since the fi-

delity term, 1 2 ∥T (n) -A (n) W ∥ 2
F , is continuous and the regularization term, equivalent to g in Theorem 2, is lsc (cf. Assumption 1 (i)). Finally, the objective function in (5.12) is bounded from below (since its values are always non-negative), and satisfies KL according to the Remark 2. Assumptions on h and g are the same as that of Theorem 1, and the same explanations provided for Theorem 1 hold also here.

The only assumption to be discussed is that of generating a bounded sequence by SFBS, which is a critical assumption in Theorem 2. [ABRS10, Remark 5] discusses about the assumptions, which guarantee the boundedness of generated sequence. For instance, the coercivity [RW09, Definition 3.25] of objective function is simply sufficient to obtain bounded sequence from

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

SFBS. A function is coercive, if it is bounded below on bounded sets and

lim inf |x|→∞ f (x) x = ∞.
It can be observed that the objective function in (5.12) is coercive.

We now review the main theoretical and practical advantages of SFBS described in Algorithm 5. Firstly, it can be applied even on non-convex and non-smooth constraints such as cardinality [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF], and this does not affect its convergence. Secondly, it does not require any critical setting of the parameters. The only parameters to be set are γ and α, where the suggested values in Algorithm 5 are almost always suitable. Thirdly, compared to stateof-the-art methods such as APG and BC-VMFB, SFBS is easy enough to understand and implement. Fourthly, contrary to BC-VMFB, SFBS works with variables in matrix form, so there is no need to vectorize loading matrices. This brings an advantage in working with large dimension tensors.

Moreover, as it is developed in the next section, SFBS can be adapted to many different constraints.

Simulation

The computer experiments described in this section can be classified in three main groups:

• Synthetic data: the decomposition of data tensor synthetically generated under the non-negative and the simplex constraint;

• Synthetic data with hidden relation8 : the decomposition of the third order moments tensor of synthetic data generated according to the generative processes described in Section 3.4.1 under the nonnegative and the simplex constraint;

• Real data: the decomposition of the third order moments tensor of a well-known text data set, namely 20 Newsgroups, which consists of

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

approximately 20000 posts on 20 topics [START_REF] Nigam | 20 Newsgroups data sets @ONLINE[END_REF], under the simplex constraint.

We compare SFBS with those described in Section 5.3, namely AO-ADMM [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF], APG [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF], FastNTF-APG [ZZZ + 16] and BC-VMFB

[VCTMM17] 9 . The software of CP decomposition via least squares, namely

Nway, from a well-known tensor Tool-box [BK + 21] is also added in our comparisons. Algorithms are either tested on artificially generated tensors, which are corrupted with additive noise, or on some artificially estimated third order moments, which contain intrinsic noise due to the lack of accuracy in estimating moments. In addition, some experiments on a real text data set (20 Newsgroups) are also added to complete this chapter.

In order to compare the performances of algorithms fairly, we require all the algorithms to iterate until either they reach a maximum predefined number of iterations (denoted by iterations max-number ) or the variation of relative objective value between two successive iterations is less than a desired small value, namely ϵ 1 . If we denote the objective value in iteration k by Ψ(k), which is the difference between the input tensor and the estimated tensor, the criterion to stop iterations is as follows:

∆Ψ(k) = Ψ(k) -Ψ(k -1) Ψ(k) ≤ ϵ 1 . (5.19)
In addition, the reported result of each experiment is averaged over several realizations of the tensors, and for each realization, all the methods are initialized by an identical set of initializations to choose the best initialization point. Therefore, it is possible that not the same initialization point will be the best for all the methods, but nonetheless, the comparison is fair, since the set of initializations is the same for all the considered methods. Therefore, the original codes of authors have been used to obtain the results reported in all figures of this section.

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

In the sequel, we mention the number of the realizations of tensors and the number of initialization points by average number and initialization number, respectively.

All computer experiments reported in this section have been executed

either on a laptop with a processor of 3.1 GHz Intel Core i5, 16 GB RAM or on a PC with a processor of 3.2 GHz Intel Core i5, 8 GB RAM, both running macOS Mojav and MATLAB 2019a.

Synthetic data

According to the considered constraints (non-negativity or the simplex set), we generate randomly (uniform distribution in the interval [0, 1]) loading matrices and a coefficient vector. Then, the noiseless tensor T o is computed via (2.6).

In order to work in a noisy context, a noise tensor, T n , with i.i.d. entries of Gaussian distribution with zero mean and unit variance, of the same size as T o , is weighted by the parameter σ and added to T o . As T n has unit variance, then the variance of σT n is σ 2 . σ is adjusted such that we reach a desired Signal to Noise Ratio (SNR) according to the following relation:

SNR = 10 log 10 1 M i,j,k T o (i, j, k) 2 1 M i,j,k σ 2 T n (i, j, k) 2 , (5.20)
where M is the total number of elements in tensors T o or T n .

Denote by T the desired tensor to be decomposed, T = T o + σT n . After decomposing T , the estimation of T o can be calculated through (2.6), as a rank-R approximation, which we call T . The relative reconstruction error is computed as follows:

ϵ( T ) = ∥ T -T o ∥ 2 F ∥T o ∥ 2 F , (5.21)
which will be reported as it is (not in the form of percentage 10 ) for all the experiments of synthetic scenarios. In addition to the reconstruction error, we will report the error of estimating loading factors but according to the points remarked in Remark 3.

Remark 3. As mentioned in Chapter 4, it is hard to assess the relative error made on loading matrices, because of the scaling and the permutation ambiguities of tensor decomposition [START_REF] Comon | Tensors: a brief introduction[END_REF]. So as to overcome these ambiguities, we report CorrIndex measure described in Chapter 4 as the performance on estimating loading factors. However, if we report CorrIndex on each loading factor separately, it would be an optimistic measure, since implicitly a specific (not common) permutation is permitted for each loading matrix.

In order to have a more reliable performance index, we report CorrIndex based on the matrix X which consists of all loading factors together. In other words, we report CorrIndex(X, X), where X, X are defined as follows

X =         A (1)
A (2) . . .

A (N )         , X =         A (1)
A (2) . . .

A (N )         . (5.22)
Nevertheless, it is worth to report also exact errors based on the methods described in 4.3.2 ( e.g. Hungarian), and compare them with the results interpreted from CorrIndex.

We should note two technical points about our practical implementations.

First, in order to compute Lipschitz constant (line 4 of Algorithm 5), one can employ [norm(W)] 2 command in MATLAB, where W is the matrix composed of constant loading matrices and is defined in line 3 of Algorithm 5. Second, as all the investigated algorithms in this chapter require the rank R, as their input, one could employ Corcondia [START_REF] Bro | A new efficient method for determining the number of components in parafac models[END_REF] to obtain an estimation of R. However, in synthetic scenarios, we are aware of real rank. In working with real text data set, some points on choosing the rank are remarked in Section 5.6.3.

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

Non-negativity constraint

Small size tensors. In all simulations of this section, ϵ 1 = 10 -20 , e = 1.9

and we repeat five times the loop of line 6 for each mode in Algorithm 5

(SFBS algorithm).

Figure 5.3 shows the relative reconstruction error of the noiseless tensor of size 10 × 10 × 10, of rank R = 6 with initialization number, average number and iterations max-number, equal to 10, 10 and 5000, respectively (cf. see page 101 for the definition of these parameters). As it can be seen at the first glance, BC-VMFB, in spite of our effort to adjust its parameters properly, results in a relative error around 1, which means 100% (remind that we always report relative error according to (5.21) which is not in the form of percentage). Although we consulted the corresponding author of BC-VMFB about the parameters of their algorithm, we found it hard to adjust several parameters of BC-VMFB. In fact, not only these parameters depend on data, but also they have a critical effect on the final result of the method. Therefore, it seems unnecessary to include BC-VMFB in the rest of comparisons in this experiment. In addition, since FastNTF-APG tries to decompose a low-rank approximation of the desired tensor, it can be expected that the performance of FastNTF-APG is worse than APG.

As a result, we remove FastNTF-APG in some comparisons of this chapter.

Although most of the methods achieve reasonable performances, APG and SFBS converge rather faster.

In Fig. 5.4, the gap between X and X via CorrIndex and Hungarian algorithms (which is an exact error) is reported for the same experiment of Fig. 5.3. As it can be seen in Fig. 5.4, the same interpretation as that of Fig. 5.3 can be concluded from CorrIndex and Hungarian. Moreover, the gap between loading factors11 and their estimations (Fig. 5.5 -Fig. 5.7) reveals the same conclusion, that is, APG and SFBS converge slightly faster than others, whereas all the methods (except FastNTF-APG and BC-VMFB) reach a reasonable performance. As the decomposition of large noisy tensors is a difficult problem to handle, FastNTF-APG is not anymore better than APG like the experiment in Fig. 5.8. In addition, unlike Fig. 5.8, Nway reaches the same level of relative error as others after convergence. Nevertheless, in this experiment, AO-ADMM and SFBS outperform others and also SFBS is slightly better.

In Fig. 5.13, the gap between X and X via CorrIndex and Hungarian algorithms (which is the exact error) is reported for the same experiment of Fig. 5.12. As it can be seen in Fig. 5.13, SFBS performs slightly better in estimating loading factors compared to the others.

Simplex constraint

In this section, we investigate practically the performance of tensor decomposition algorithms under the simplex set constraint over all the columns of all loading factors and over λ. In the simulations of this section, we set ϵ 1 = 10 -20 and repeat five times the loop of line 6 for each mode in Algorithm 5 (SFBS algorithm).

In this section, we compare SFBS with AO-ADMM, since, firstly, the include the information of the estimation of λ, we compare the performance of algorithms in the estimation of λ in Fig. 5.17, which reveals the same conclusion as the other plots of this experiment.

Noisy cases. The same conclusion can be drawn as before, that is, SFBS performs and converges better and faster than AO-ADMM. It can be seen that the relative error in Fig. 5.20 increases slightly, but note that decreasing monotonically is essential for the relative reconstruction error (or objective function) not in the estimation of loading factors or λ, which is the case in all the experiments such as the one expressed in Fig. 5.18.

Synthetic data with hidden relations

In order to follow much more easily the computer experiments of this section, it is better to review briefly the notation of the data mining problem described in Section 3.3:

• h: a hidden variable (topic) encoded into a discrete variable taking K possible integer values, say in H = {1, 2, . . . , k, . . . , K} with the probability φ(k) = Prob(h = k)

• Ω = {u 1 , . . . , u d , . . . , u D }: the encoded dictionary of cardinality D for the multi-view variables (words), x

• f k (d) = Prob(x = u d |h = k): the conditional probability of each word u d of dictionary Ω, given a particular topic, h = k • a k : a vector of dimension D contains the values of f k (d) for all d • A = [a 1 , . . . , a k , . . . , a K ].
The goal is to estimate the probabilities A and φ by means of the decomposition of the third order moments tensor (T ) of the observed multi-view variables (words), which can be employed for some targeted data mining tasks such as unsupervised clustering. As it is described in Section 5.4, some algorithms such as Robust tensor power method require also the second order moments (P ) for the decomposition of T . In addition, some methods of moment estimation, such as simple and standard averaging, have been expressed in Section 3.4.

The drawbacks of Robust tensor power method [AGH + 14]

In this section, we show experimentally the drawbacks of Robust tensor power method, namely the negative estimated probability (which is also reported in our paper [START_REF] Sobhani | Data mining with tensor decompositions[END_REF] but we shall explain it completely in the sequel)

and the sensitivity to the additive noise. Figures 5.21 power method occasionally are not acceptable, since they include negative values. Table 5.2, Fig. 5.23 and Fig. 5.24 experimentally prove that Robust tensor power method is very sensitive to the additive noise.

Negative estimated probabilities. In the simulation of Fig. 5.21 and Middle and right column in Fig. 5.21 correspond to the estimation of the columns of matrix A with the non-negative AO-ADMM and Robust tensor power method, respectively. In order to compare the performances, the plots of each row of Fig. 5.21 should be compared with each other. In all the rows of Fig. 5.21 except the first row, the negative values can be seen in the estimation of Robust tensor power method, which are not acceptable as the estimated probabilities for the targeted data mining task.

The performances in estimating vector φ have been compared in Fig. 5.22, where the figures in the middle and at the bottom correspond to the nonnegative AO-ADMM and Robust tensor power method, respectively. It is clear that the result of the non-negative AO-ADMM is closer to the original DECOMPOSITIONS vector φ.

Sensitivity to additive noise. In the sequel, we shall experimentally show that Robust tensor power method does not perform well in presence of additive noise. In order to assess this, we generated some arbitrary matrices A 10×3 and vectors φ 3×1 . Then, we calculated their original corresponding third and second order moments, i.e. T and P (cf. (3.5) and (3.6)). At the end, we applied several tensor decomposition algorithms, including Robust tensor power method, on both the noiseless (Table 5.2) and the noisy (Fig. 5.23 and Fig. 5.24) versions of T and P .

Table 5.2 reports the relative reconstruction error and the relative exact error in estimating A and φ (measured with Hungarian algorithm) from noiseless T and P . These results averaged over 200 realizations of A and φ14 . The considered algorithms in this experiment are as follows:

1 -Power method: Robust tensor power method 2 -Projected Power method: as the result of Robust tensor power method is not in the simplex set, we modified this method by projecting its result to the simplex set.

3 -Simplex AO-ADMM: AO-ADMM with the simplex set constraint 4 -Simplex SFBS: SFBS with the simplex set constraint 5 -Symmetric Simplex SFBS: unlike Simplex SFBS, we take into account the symmetric property of the third order moments tensor. To be more precise, in calculating the objective function for the criterion to stop the iterations of SFBS, we utilized an average of all estimated loading factors. For the last three methods (i.e. Simplex AO-ADMM, Simplex SFBS and Symmetric Simplex SFBS), we set initialization number to 20, iterations max-number to 1000, ϵ 1 = 10 -20 and e = 1.9. In addition, in order to have a fair comparison with Power method, we set "M 1 = initialization number = 20" and "M 2 = 5 × iterations max-number" in Algorithm 4, where 5 is due to the five times of doing the loop of line 6 in Algorithm 5 (the same loop in AO-ADMM is also done five times).

It can be inferred from Table 5.2 that all the considered methods perform well in noiseless case, specially Projected Power method. These results show that our implementation of Power method is reliable.

Although Power method and Projected power method perform very well in the noiseless case, Fig. 5.23 and Fig. 5.24 show that these algorithms are very sensitive to additive noise in T and P . Note that in these figures SFBS, SFBS-Symm, ADMM, Power and Power-proj refer to Simplex SFBS, Symmetric Simplex SFBS, Simplex AO-ADMM, Power method and Projected Power method, respectively.

All the settings in the experiment of Fig. 5.23 and Fig. 5.24 are the same as that of Table 5.2, except that T and P are perturbed by additional Gaussian noise at a specific SNR (cf. (5.20)). expresses that even by projecting the result of Power method to the simplex set (Projected Power method), the performance of Power method cannot become acceptable.

The relative error of estimation of A and φ are depicted in Fig. 5.24 (Left) and Fig. 5.24 (Right), respectively. As it is demonstrated in Fig. 5.24, Power method and Projected Power method does not perform as well as others such as SFBS. Moreover, Simplex SFBS and Symmetric Simplex SFBS perform a bit better than AO-ADMM. Note that unlike reconstruction error, the relative error of the estimation of A and φ is limited in the range of [0, 2], since we normalize the φ and the columns of A with L 2 norm before measuring the performance index. To carry out this sample estimate, some moment estimators such as simple averaging are reviewed in Section 3.4.3, as well as our proposed estimator, standard averaging, in Section 3.4.4.

Hence, it is essential to check the consistency (cf. see Section 3.4.2 for a definition of consistency) of generated data along with moment estimator in order to make sure that sample moments indeed converge to the true joint probabilities.

The effect of constraint on the performance. In this experiment, we experimentally show the effect of considering a constraint on the precision of probability estimation. In other words, we shall compare the performance of non-constrained tensor decomposition such as Robust tensor power method with the constrained one, such as AO-ADMM and SFBS, in terms of their Each black circle (resp. red cross) corresponds to the discrepancy between sample third order moments (resp. sample second order moments) and true moments for a specific corpus of documents with a particular corpus size, N c . The median among corpuses of documents is also plotted in solid line for every corpus size.

Data generated by the generative process described in 3.4.1.1 are consistent for N c > 10000. Sample moment computation is carried out via simple averaging. The relative error is reported according to (5.21), and it is not in the form of percentage.

choice of A and φ):

A =                    0.
                   , φ =        0.256 0.163 0.201 0.380       
.

In order to compare the performance of both methods, we run this simulation for increasing the corpus size (note that the larger the corpus size, the fewer the number of documents, to keep a computational load constant).

To be more precise, by keeping the values of parameters D = 8, K = 4, we increase exponentially the value of N c from 2 14 to 2 17 . For each value of N c , we generated several synthetic data sets of size N c in the way that has been explained in Subsection 3.4.1, and then the methods have been applied to the obtained data to compare relative errors as a function of corpus size.

The results are reported in Fig. 5.26; each black circle (resp. red cross) corresponds to the error of Robust tensor power method (resp. non-negative AO-ADMM) in estimating φ or A for a specific corpus of documents with a particular corpus size, N c . To ease comparison, the median (resp. standard deviation) among corpuses of documents is also plotted in solid (resp. dotted) line for every corpus size.

As Fig. 5.26 shows, the error of constrained method such as non-negative AO-ADMM converges faster to zero as the corpus size increases, and unlike Robust tensor power method, non-negative AO-ADMM appears to be much more robust because of its smaller standard deviation of the relative error. In Table 5.3 the results obtained for a particular value of K and a few values of D, averaged over 20 independent trials. These extensive computer experiments confirm the robustness and superiority of constrained algorithms such as non-negative AO-ADMM compared to Robust tensor power method.

Improve performance by standard averaging and simplex set constraint. We can modify the experiment of Fig. 5.26 from different aspects.

As expressed in Section 3.4.5, the sample moment estimated by standard averaging requires a smaller corpus size compared to simple averaging to converge to the true moments. In addition, as we seek the probabilities, which belong to the simplex set, the decomposition would be more precise if it is carried out under the simplex set constraint for φ and the columns of

A.
Therefore, one can feed decomposition algorithms mentioned in page 118 by the sample moments approximated via standard averaging. Moreover, it is preferred to perform constrained algorithms such as AO-ADMM and SFBS under the simplex set constraint. These modifications have been done in the experiment of Fig. 5.27.

The experiment of Fig. 5.27 is performed for K = 3 (hidden variables (topics)), D = 10 (multi-view variables (words)) and various corpus size up to N c = 2 10 . The synthetic data for this experiment are generated via the generative process described in Subsection 3.4.1.2, which is based on multinomial distribution. The considered corpuses are consist of documents with minimum and maximum length of 3 and 100, respectively. We set initialization number to 20, iterations max-number to 1000, ϵ 1 = 10 -20 and e = 1.9.

In addition, in order to have a fair comparison with Power method, we set "M 1 = initialization number = 20" and "M 2 = 5 × iterations max-number" in Algorithm 4, where 5 is due to the five times of doing the loop of line 6 in Algorithm 5 (the same loop in AO-ADMM is also executed five times). Let us remark two other points that can be inferred from Fig. 5.27.

First, by comparing Fig. 5.26 and Fig. 5.27, we can observe easily the effect of moment estimators. In Fig. 5.27, which is based on standard averaging, a corpus of size 1024 is sufficient for the moment consistency, hence, for an acceptable estimation of probabilities. Nevertheless, in Fig. 5.26, a corpus of size N c > 10000 is necessary to have a reasonable performance, since the moment estimator is simple averaging. Second, for the smallest considered corpus size in Fig. 5.27, i.e. N c = 32, the SNR = 7.62, which shows that our assumption on the range of SNR in the synthetic data with hidden relation (Fig. 5.23 and Fig. 5.24) was reasonable.

Real data

In this section, we describe our experiment on a part of a well-known text data set, namely 20 Newsgroups, which consists of 11314 posts on 20 topics available online [Lan95a, Lan95b, Nig00].

It is vey common to do some pre-processing steps in order to extract the keywords of the corpus described above. standard steps of pre-processings are as follows: removing stop words (including me, you, a, the, which, where, ...), tokenizing (transform a text into a list of words), removing punctuations and unnecessary characters, lemmatizing (replacing with the infinitive, e.g. replacing "took" with "take"). At the end, according to the frequency of each word (called term-document frequency) in the corpus, a portion of words are kept as the representatives of the corpus also known as vocabulary or dictionary.

By doing above mentioned pre-processing15 steps on the whole data set and by keeping the words with more than 20% of the term-document frequency, we are left with a dictionary of size 14 × 14 (D = 14 words).

Since, Power method and its variants cannot handle over-complete cases (K > D), we are forced to choose K = 14 topics rather than K = 20, though K = 20 is more suitable according to 20 topics exist in the considered data set (N c = 11314 posts on 20 topics). Therefore, the resulted bag-of-words (cf. (3.9)) is a matrix of dimensions 14 by 11314.

After obtaining the bag-of-words, we fed the tensor decomposition algorithms listed in page 118 with the estimated P and T via standard averaging.

In addition, we set the input rank of all the algorithms to K = 14 due to the limitation of Power method and its variants. Eventually, by each algorithm, a pair of estimated probabilities (φ 14 , A 14×14 ) is acquired.

Although in this corpus, the number of documents for all the topics are almost the same (which means that φ is expected to have a distribution close to the uniform distribution), the estimated φ by Robust power method, i.e. φ Power , is not only non-uniform, but also it is not a probability distribution at all, since it does not lie in the probability simplex: As it can be observed in φ Power , its tenth element dominates the others which could be the consequence of rounding errors. Even if one projects φ Power to the simplex set as mentioned before, and obtains φ Power-proj , it would be meaningless according to the fact that "φ is expected to have a distribution close to the uniform distribution", since we have:

φ Power = [0.
φ Power-proj = 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 T .

On the other hand, the estimated φ by SFBS, i.e. φ SFBS , is as follows:

φ SFBS = [0.10, 0.08, 0.09, 0.09, 0.10, 0.04, 0.02, 0.05, 0.06, 0.01, 0.06, 0.09, 0.08, 0.12] T , which is closer to uniform distribution than φ Power-proj . This again proves that applying constrained tensor decomposition is more reliable than nonconstrained one.

In Section 5.6.2.1, we investigated two drawbacks of Robust tensor power method, namely negative estimated probabilities and sensitivity to the additive noise. Nevertheless, in this experiment, we detect a much more important drawback in employing Power method for real text data set, that is, the dominant estimated element.

Although Projected power method resolves the problem of negative estimated probabilities by projecting the estimated probabilities to the simplex set, it is incapable to solve the problem of the dominant estimated element.

Therefore, utilizing non-constrained tensor decompositions such as Power method and its variants is not a good solution in practice, however these kinds of methods are considered as fast solutions for estimating probabilities.

In order to be able to compare the performances of tensor decomposition algorithms on this real text data set, we need the ground truth values of the probabilities of topics, i.e. φ, and the conditional probabilities of each word given the topic, i.e. A. To do this, for the sake of convenience, we worked with a portion of data set introduced at the beginning of this section. We selected documents about four different topics, namely "computer graphics" (labeled by comp.graphics), "baseball" (labeled by rec.sport.baseball), "cryptography" (labeled by sci.crypt) and "Christianity" (labeled by soc.religion.christian). Then, by applying above mentioned pre-processings on this portion of data, which contains 2375 documents, and by keeping the words with term-frequency between 20% and 50%, we had a dictionary with D = 17 words. As we desired to employ DECOMPOSITIONS Ruffini's estimator in this experiment, it was necessary to eliminate the documents including less than three words. By doing that, 1690 documents remained in the corpus. So as to calculate the ground truth values of φ, we counted the number of documents belong to each topic (i.e. "computer graphics", "baseball", "cryptography" and "Christianity"). Then, the ground truth of φ is simply the division of the number of documents in each topic by 1690. To calculate the ground truth of A, for each word in the obtained dictionary, we counted its occurrences in the documents belong to a particular topic, and then we divided it by the total number of words used in those documents (which belong to that particular topic).

As in the previous experiment, after obtaining the bag-of-words, we fed the tensor decomposition algorithms listed in page 118 with the estimated P and T via the Ruffini's estimator. We set the input rank of all the algorithms to K = 4, and eventually, we computed the performance of each algorithm by comparing estimated probabilities ( φ 4 , A 17×4 ) with the ground truth values of φ and A.

The performances of each algorithm of this experiment in estimating probabilities φ and A based on CorrIndex are reported in Table 5.4. As it can be seen in Table 5.4, the performances of Power method and Projected Power method in estimating φ are much worse than other constrained algorithms.

However, their performances in estimating A are a bit better than others.

In addition, all the constrained algorithms such SFBS perform properly in both estimating φ and A.

The fact that Power method performs better in estimating A than in estimating φ is probably due to the different manner of updating φ from updating A in Power method. As it can be observed in line 6 of Algorithm 4, the L 2 norm of the columns of A are stored as the estimation of φ. In this way, a normalized matrix is reported as an estimation of A, which is protected relatively from rounding errors. However, the non-normalized φ may suffer from rounding errors. As the words in the obtained dictionary of this experiment 16 are not informative enough, it seems impossible to do an unsupervised clustering based on the estimated probabilities A and φ. However, extracting informative words depends on employing proper pre-processing steps, which is out of the scope of this thesis.

Discussion

According to the performed experiments and theoretical arguments, in this section, we wrap up the advantages and disadvantages of each method compared to our proposed method, SFBS.

Concerning constrained algorithms, a theoretical drawback of AO-ADMM is that the convergence of AO is not guaranteed because of ignoring the proximal regularization (cf. Subsection 5.3.2). In addition, for the non-convex constraint ℓ 0 , the convergence of ADMM has not yet been proved, whereas the complete convergence of SFBS is studied in Section 5.5.4. However, its 16 also, come, even, find, get, give, go, good, make, need, nntp-poste, people, see, take, time, way, well.

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

performance is as well as that of SFBS in most cases, such as noiseless and noisy scenarios.

By comparing APG and SFBS, we can conclude that SFBS performs better than APG in noisy cases (cf. Fig. 5.9), while it is also computationally less expensive, due to not calculating any coefficient, unlike APG, which is based on FISTA. Moreover, some methods such as BC-VMFB forces the user to adjust several parameters that are very effective on the result, while the only parameter of SFBS is e, for which we provided a proper range ([1.5, 1.9]).

Furthermore, unlike other constrained algorithms, we discussed how SFBS can handle a variety of constraints, namely non-negativity, the simplex set (for all loading factors plus the vector of coefficient), sparsity (ℓ 0 pseudonorm and ℓ 1 norm).

Concerning non-constrained algorithms, SFBS is much more robust against additive noise than Nway or Power method. Moreover, Power method and its variants (such as SS-HOPM) does not have any convergence guarantee for an odd-order tensor (including the tensor of third order moments), whereas the complete convergence of SFBS is studied in Section 5.5.4.

Focusing on data (text) mining based on the method of moments and by means of tensor decomposition, there is an intrinsic noise due to the discrepancy between true moments and sample moments. Therefore, Power method and its variants, which are not robust enough to the noise will face a problem. Moreover, it is also possible that Power method or its variants return non-acceptable negative values for estimated probabilities. Even if one projects the result of Power method to the simplex set (i.e. Projected power method), its performance in estimating φ would be the worst among others (cf. Fig. 5.27 (bottom)). More importantly, as it is shown experimentally in Section 5.6.3, another critical drawback of Power method in face of real text data set is the dominant estimated element (i.e. a very large element in φ, which drops out φ from the simplex set), and this problem cannot be solved with projection to the simplex set by Projected Power method.

In addition, in order to decompose the tensor of third order moments,

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

Power method requires also second order moments for the whitening step.

Furthermore, since after whitening, Power method performs a kind of orthogonal tensor decomposition, hence, one of its limitations is that it can only handle under-complete cases (i.e. when the number of hidden variables, K, is smaller than that of multi-view ones, D).

However, SFBS is robust against additive noise (especially in estimating φ, cf. Fig. 5.27 (bottom)), returns the estimated probabilities in the simplex set thanks to its constraint, requires only the tensor of third order moments as input, and is able to handle over-complete cases (K > D), which is necessary in some applications mentioned in [START_REF] Anandkumar | When are overcomplete topic models identifiable? uniqueness of tensor Tucker decompositions with structured sparsity[END_REF]. The only drawback of constrained algorithms such as SFBS compared to the Power method (and especially SVTD, which is faster than the Power method) is that they are more time consuming than Power method.

As a last point, the role of moment estimator should always be taken into account, since it has a significant effect on decreasing the intrinsic noise and decreasing the required data for consistency. Even the best decomposition algorithm may not return acceptable results without being provided a good estimation of moments as its input.

Conclusion and perspective

In this chapter, we investigated constrained and unconstrained tensor decompositions. As mentioned before, the main goal is to estimate probabilities by decomposing the third order moments tensor, and in the literature, usually some unconstrained tensor decompositions such as Robust tensor power method are employed for this purpose.

We showed theoretically and experimentally that such a unconstrained algorithms are not proper, since they probably return negative values, which are not acceptable as probabilities (in this vein, we proposed Projected Power method to solve this problem by projecting the results of Power method to the simplex set); or they are very sensitive to additive noise and are unable

CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

to handle over-complete cases, where the number of hidden variables is more than the number of multi-view ones. More importantly, we experimentally showed that in employing Robust tensor power method on real text data set, a very challenging problem arises, namely the dominant estimated element drops out estimated probability vector (φ) from the simplex set. Moreover, this issue cannot be solved by projecting onto the simplex set.

On the other hand, we showed that constrained tensor decompositions return more acceptable estimations, since some constraints are taken into account, such as non-negativity and simplex set. We also introduced our constrained tensor decomposition, SFBS, which is based on Forward-Backward Splitting.

We compared theoretically and experimentally SFBS with other algorithms, and we concluded that SFBS performs either better than or as well as AO-ADMM. Furthermore, there exists a complete convergence proof for SFBS, which is not the case for AO-ADMM due to ignoring the proximal regularization (which is essential for the convergence guarantee of AO) and due to the lack of theoretical proof for the non-convex constraint ℓ 0 .

We also explained that SFBS performs better than APG in noisy cases and computationally costs less than APG because of not calculating any coefficient, unlike APG. In addition, we discussed that not only BC-VMFB performs poorly in decomposing tensors, but also it forces the user to adjust several parameters, which are so effective on the result of this method. Yet, SFBS consists of just one parameter (e), which is easy to set according to the convergence condition and the spectral norm of the corresponding matrix (e.g for the experiments of this section we proposed a range of [1.5, 1.9]).

Compared to unconstrained algorithms such as Robust tensor power method, not only SFBS performs better, but also it has the theoretical convergence guarantee. Moreover, concerning the simplex set constraint, one can assure that probabilities estimated by SFBS always lie in the probability simplex. Further, SFBS is much more robust against additive noise, which is a very important feature in data (text) mining where there exists CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS intrinsic noise due to the discrepancy between true moments and sample moments. However, non-constrained algorithms, especially SVTD, converge much more faster than iterative algorithms such as SFBS.

In addition to the intrinsic noise, working with large dimension third order moments tensors is another big challenge of data mining by means of tensor decomposition. Recall that the dimension of this tensor depends on the number of multi-view variables (words) in the dictionary, which can be very large in a real world text mining problem (around 5000). Unfortunately, most of iterative constrained algorithms such as SFBS are very time consuming, hence, intractable in practice.

As the most time consuming part in tensor decomposition algorithms is the calculation of Khatri-Rao product of loading matrices, several methods have been proposed to reduce this computational complexity by sketching [START_REF] Keriven | Sketching for large-scale learning of mixture BIBLIOGRAPHY models[END_REF] the data required for the calculation of Khatri-Rao product in each iteration. One of the most successful methods in this vein is CPRAND [START_REF] Battaglino | A practical randomized cp tensor decomposition[END_REF], which is a non-constrained CP decomposition and is proper for large dimensions, since it employs sketching data for the calculation of Khatri-Rao product, called sample Khatri-Rao.

Recently, a constrained CP decomposition for large dimensions based on Forward-Backward Splitting and the sketching concept has been introduced in [FIW + 20]. Although, it is claimed that this method is capable of handling various types of constraints, but we find it inefficient for the simplex set constraint, in spite of using the implementation provided by the authors.

Therefore, a future extension of SFBS would be employing CPRAND idea (sample Khatri-Rao) in each iteration of SFBS. In this way, we will have a constrained CP decomposition dedicated for large dimension tensors, which handles properly constraints such as the simplex set constraint.

Conclusion and Perspectives

Vectors and matrices as one-way and two-way arrays, respectively, can be generalized to multi-way arrays or tensors. Tensors are known to be useful to identify parameters, thanks to the mild uniqueness conditions required in Canonical Polyadic tensor decomposition (3.6).

In this thesis, we focused on a particular data model, called single topic (cf. Fig. 3.1), consisting of some multi-view variables which are related to each other via a hidden/latent variable. A document can be described properly by this model when its words role as some multi-view variables, which are related to each other according to the topic of that document as a hidden variable.

Estimating the probabilities of hidden variables, and also the conditional probabilities of multi-view ones, are in the center of interest, since they can be utilized for an ultimate data/text mining task such as unsupervised clustering. On the other hand, these probabilities correspond to the tensor of third order moments of observed data (e.g. words) via a particular tensor format, called Canonical Polyadic (CP); see (3.6) for details.

In this thesis, generally, we investigated CP decomposition applications in estimating the probabilities of hidden variables and the conditional probabilities of multi-view variables.

In chapter 3, we explained that there are two main steps in estimating the probabilities via a tensor approach: firstly, estimate the third order moments tensor by observing data which obeys the single topic model, secondly, decompose the estimated third order moments tensor by means of proper tensor decomposition algorithms to obtain some reliable and accu-rate estimations of probabilities that can be used efficiently for an ultimate data/text mining task such as unsupervised clustering.

At the end of chapter 3, we reviewed some existing moment estimators, and criticized simple averaging, as a relevant moment estimator in the literature, by experimentally showing its weakness in providing moment consistency even on sufficiently large size corpuses.

As our moment estimator, we proposed standard averaging, which performs much better than simple averaging, due to taking into account all the words of each document as well as the length of each document. In this way, we control the importance of each document in the estimation. Further, compared to the state-of-the-art, standard averaging is based on a simpler concept (weighted averaging), and is computationally less expensive in terms of number of multiplications for estimating third order moments. In addition, the number of required additions for estimating second order moments does not depend on the dictionary size.

In Chapter 4 and 5, we discussed about the second stage, i.e. tensor decomposition algorithms, and also about required performance indices to evaluate these different methods of decompositions. Recall that permutation and scale ambiguities are inherent in tensor representations, by definition of tensors [START_REF] Comon | Tensors: a brief introduction[END_REF]. Therefore, in order to cancel the effect of these ambiguities in measuring the performance, one needs to employ proper performance indices.

Chapter 4 is devoted to performance indices in which we carried out a critical survey about the state-of-the-art. We showed theoretically and experimentally that existing indices are either greedy (optimistic or pessimistic) or computationally expensive. Moreover, we proposed a new performance index belonging to the class of invariant indices, called CorrIndex, whose upper and lower bounds are easy to interpret, while being computationally cheap.

We compared tensor decomposition algorithms in Chapter 5 by means of CorrIndex.

In Chapter 5, we investigated constrained and unconstrained tensor de-compositions to discover their advantages and disadvantages especially in decomposing the tensor of third order moments to estimate probabilities.

We experimentally showed that although Power method and its variants are relevant methods for this purpose, they are sensitive to additive noise and not capable of handling over-complete cases, where the number of hidden variables is larger than multi-view variables. More importantly, the probabilities estimated by these unconstrained algorithms may include negative values or dominant elements, which exclude estimated probabilities out of the simplex set. One may think of projecting the estimated probabilities to eliminate negative values but this solution is not useful for the problem of dominant estimated element.

On the other hand, we proposed to employ constrained tensor decompositions such as AO-ADMM under non-negative or simplex set constraints in order to have much more accurate estimations of probabilities. In this vein, we proposed our constrained tensor decomposition, called Simple Forward-Backward Splitting (SFBS), using Forward-Backward Splitting, a minimization procedure based on proximal concept. The term "simple" shows two aspects of our algorithm: first, compared to state-of-the-art, it is less computationally costly, second, its algorithm is much easier to understand and implement. Furthermore, we described how SFBS can handle a variety of constraints. Moreover, a complete convergence analysis is provided, which is not possible for AO-ADMM. More importantly, SFBS performs better than state-of-the-art in noisy scenarios while being computationally cheap.

In summary, we investigated each required step of probability estimation via a tensor approach. We studied theoretically and experimentally the challenges of state-of-the-art in moment estimations and tensor decompositions. Moreover, we proposed our moment estimator (standard averaging), our performance index (CorrIndex) and our tensor decomposition algorithm (SFBS), which bring some advantages over existing methods and indices.

As mentioned at the end of Chapter 5, the main disadvantage of constrained tensor decomposition algorithms for large dimension (typically of the order D > 50) is that they are time consuming due to the calculation of Khatri-Rao product of loading matrices. Several methods have been proposed to reduce this computational complexity by sketching the data required for the calculation of Khatri-Rao product in each iteration. One of the most successful methods in this vein is CPRAND [START_REF] Battaglino | A practical randomized cp tensor decomposition[END_REF], which is a non-constrained CP decomposition and is proper for large dimensions, since it employs sketching data for the calculation of Khatri-Rao product, called sample Khatri-Rao.

Therefore, a future extension of SFBS would be to follow the CPRAND idea (sample Khatri-Rao) in each iteration of SFBS. In this way, we would have a constrained CP decomposition dedicated to tensors of large dimension, which handles properly constraints such as the simplex set constraint.

It is worth mentioning that recently some constrained decomposition algorithms for large dimension tensors are proposed in [FIW + 20]. However, we found them improper under the simplex set constraint. Hence algorithms in [FIW + 20] could be adapted for this constraint, which is important in estimating probabilities.

The procedure described in Sections 5.5.2.3 and 5.5.2.4 is useful when sparsity constraints are applied on each element of loading factors. However, if sparsity constraints are required for each column of loading factors (as in co-clustering application [START_REF] Evangelos | Coclustering as multilinear decomposition with sparse latent factors[END_REF]), ℓ 0 and ℓ 1 should be replaced by mixed norms [START_REF] Kowalski | Sparse regression using mixed norms[END_REF] to ensure the sparsity of each column.

Although we focused on single topic model in this thesis, the content can be extended to other data models such as LDA [START_REF] David M Blei | Latent dirichlet allocation[END_REF] and Markov models [AGH + 14].

A Lipschitz constant of the gradient of the fidelity term in (5.12)

Let us define h(A (n) ) and g(A (n) ) as follows:

h(A (n) ) ≜ 1 2 ∥T (n) -A (n) W ∥ 2 F g(A (n) ) ≜ i C A (n) (A (n) ),
where W = (A (N ) ⊙ . . . ⊙ A (n+1) ⊙ A (n-1) ⊙ . . . ⊙ A (1) ) T . Note that the objective function of (5.12) can be written as h(A (n) ) + g(A (n) ).

Since h(A (n) ) has a quadratic form, the "Lipschitz" constant of its gradient can be calculated. The gradient of h(A (n) ) is computed as follows:

h(A (n) ) = 1 2 ∥T (n) -A (n) W ∥ 2 F = 1 2 trace{(T (n) -A (n) W ) T (T (n) -A (n) W )} ⇒ ∇ A (n) h(A (n) ) = -W T (n) T + A (n) W W T
In calculating the gradient, the following relations have been used [START_REF] Duchi | Properties of the trace and matrix derivatives[END_REF]:

∇ A trace{AB} = B T , ∇ A trace{ABA T C} = C T AB T + CAB.
Now, the Lipschitz constant of ∇ A (n) h(A (n) ) can be calculated as follows:

B Résumé en Francais Sauvegarder l'index des mots du document "n".

B.2.3 Résultats de la simulation

Toutes les expériences informatiques rapportées dans cette section ont été exécutées soit sur un ordinateur portable avec un processeur Intel Core i5 de 3,1 GHz, 16 Go de RAM, soit sur un PC avec un processeur Intel Core i5 de 3,2 GHz, 8 Go de RAM, tous deux exécutant macOS Mojav et MATLAB 2019a.

B.2.3.1 L'effet de la taille du corpus sur la cohérence du moment

Comme mentionné précédemment, le simple calcul de la moyenne ne permet pas d'obtenir une cohérence des moments acceptable et une façon naïve de compenser cet inconvénient est d'augmenter la taille du corpus (le nombre de et en substituant les valeurs de |a H i âj | de C ij , on obtient 3-0.85-0.9-0.7 = 0.55, ce qui est inférieur à l'erreur exacte, 3 -0.8 -0.9 -0.7 = 0.60 (l'erreur exacte est donnée dans la section 4.3.3 avec la permutation optimale). Cet exemple montre que cet algorithme produit une matrice P p qui peut ne pas être une permutation. Cet indice est toujours optimiste (cf. Section 4.3.1.1).

Approche gloutonne de [CLDA09, CKAC14] Afin d'éviter une affectation non acceptable, après avoir détecté la valeur maximale de chaque colonne de C, sa ligne et sa colonne peuvent être supprimées pour le reste de l'algorithme. Il s'agit d'une approche gloutonne, puisque l'indice dépend de l'ordre de choix des valeurs maximales. Par exemple, si cet algorithme glouton est appliqué à la matrice C exprimée en (B.12), l'affectation résultante sera (â 1 , a 2 ), (â 2 , a 1 ), (â 3 , a 3 ) à condition que les colonnes soient balayées de gauche à droite. En revanche, si les colonnes sont balayées dans le sens inverse, l'affectation sera (â 1 , a 1 ), (â 2 , a 2 ), (â 3 , a 3 ). Par rapport à l'indice optimiste, l'erreur produite par cette approche gloutonne en balayant de gauche à droite, 3 -0.85 -0.3 -0.7 = 1.15, est plus grande que l'erreur exacte, 3 -0.8 -0.9 -0.7 = 0.60, tandis qu'en balayant de droite à gauche, l'erreur rapportée est égale à l'erreur exacte 0.6. Par conséquent, l'erreur signalée est toujours supérieure ou égale à l'erreur exacte sur la base de l'affectation optimale; cet indice est donc pessimiste. est défini comme l'opérateur de proximité de la fonction f ∈ Γ 0 (R N ), et il est désigné par prox f (x).

Dans de nombreuses applications de traitement du signal, la fonction de coût à minimiser est la somme de deux fonctions dont l'une est généralement non-différentiable ou même non-convexe. En suivant une approche proximale, ces types de problèmes peuvent être résolus à l'aide d'un algorithme particulier appelé Séparation Avant-Arrière.

y k = z k -β -1 ∇h(z k ) 4:
x k+1 = prox β -1 g (y k ) 5:

t k+1 = 1 + 4t 2 k + 1 2 6: λ k = 1 + t k -1 t k+1
7: for g = 0, 1, 2, . . . do 7:

z k+1 = x k + λ k (x k+1 -x k
Y = A (n) g -γ(A (n) g (W W T ) -T (n) W T ) 8: A (n) g+1 = A (n) g + α g (proj C A (n) (Y ) -A (n) g ) 9:
end for 10:

end for 11: until Un certain critère de fin l'équation à point fixe suivante :

A (n) = prox γi C A (n) {A (n) -γ(A (n) W W T -W T (n) T )}.

viii

  Mots clés-Décomposition tensorielle, Fouille de texte, Clustering non-supervisé, Variable cachée/latente, Moment de troisième ordre, Fractionnement avant-arrière, Opérateur proximal, Indice de performance, Ambiguïté d'échelle et de permutation Contents The tensor A of dimension 2 × 2 × 2. . . . . . . . . . . . . . . 3.1 The mixture of unigrams or single-topic model with its multiview variables (x i ) and their corresponding hidden variable (h). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  figure, le nombre de corpus utilisés est fonction de leur taille: plus la taille du corpus est grande, moins il y a de corpus; l'intérêt de ce choix délibéré est de maintenir constante la longueur totale du bloc de données (et donc la charge de calcul) pour chaque point en abscisse. . . . . . . . . . . . . . . . 154

5. 1

 1 Proximity operator of functions used in this thesis . . . . . .5.2 Generating some arbitrary matricesA 10×3 and vectors φ 3×1 , then calculating their corresponding T and P , averaging the results of decomposing noiseless T over 200 realizations of A and φ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxv xxvi LIST OF TABLES 5.3 The median and standard deviation of error in estimating A and φ, with a corpus of size N c = 2 17 ≈ 1.2 × 10 6 , and K = 4.

  . . . . . . . . . . . . . . 135 B.1 Comparaison numérique des méthodes de mesure de la distance entre A 150×100 avec une cohérence mutuelle γ = 0, 75 et sa version bruitée permutée A avec SNR = -1, 76 dB en moyenne sur 50 de réalisations. L'indice dans les cinq premières lignes du tableau est l'erreur relative. En revanche, les quatre derniers indices du tableau sont définis différemment et ne sont donc pas comparables. . . . . . . . . . . . . . . . . 163

0

 0 product ∥.∥ p norm p ∇f The gradient of a differentiable function ∂f The sub-gradient of a function f β The Lipschitz constant of the gradient of the differentiable function h in Forward-Backward Splitting theorem Γ The class of lower semi-continuous functions γ e mapping from L to Ω δ Dirac delta function ϵ 0 The gap between A and AP σ Λ Probability of topic k φ K The vector of probabilities of dimension K: [φ(1), φ(2), . . . , φ(K)] T φ The estimation of φ xxxi Φ(k)

ηηG

  estimation of A b n The bag-of-words of document n C Correlation matrix d word index in a dictionary D Number of multi-view variables (words) in a dictionary dom(f ) The domain of a function f E Expectationf k (d) Conditional probability of word d when the topic is k F k (d)The cumulative distribution of words:d i=1 f k (i) F D×KThe matrix of cumulative distribution of wordsf hThe hard-thresholding function with a threshold η f sThe soft-thresholding function with a threshold η Core The number of words in document nN cThe size of a corpus (a set of texts)N Tensor order prox f (x)The proximity operator of a function f at point x proj S (x)The projection of a vector x onto a closed convex set S

  th columns of N loading matrices (factors) A (1) , A (2) , . . . , A (N ) , respectively.The N matrices A (1) , A (2) , . . . , A (N ) are called mode-1, mode-2, . . . , mode-N loading matrices, respectively, since their columns are responsible for the construction of the first, second, . . . , N th dimension of T [CMDL + 15].

  probability of each word u d of the dictionary Ω, given a particular topic, h = k • a k : a vector of dimension D containing the values of f k (d) for all d • A = [a 1 , . . . , a k , . . . , a K ].

CHAPTER 3 .

 3 DATA MINING AND TENSORS 3.4.1.1 Generative process based on the cumulative distribution To generate a data set, we need first to define the distributions φ(k) and f k (d) for all (k, u d ) ∈ H × Ω. The values of φ(k) are stored in a K-dimensional vector φ. Similarly, the values of f k (d) are stored in a D × K matrix A. But it is actually useful to use their cumulative distributions, Φ(k) = Prob(h ≤ k) and F k (d) with an appropriate encoding, as we shall see. The values of Φ(k) are stored in a K × 1 vector Φ obtained from φ by cumulated sum of its entries. Similarly, the values of F k (d) are obtained by cumulated sum of the entries of matrix A over first dimension and stored in a D × K matrix F . Our generative algorithm goes along the following lines:

( 3 .

 3 10) to (3.14) express that in spite of simple averaging, Zou's estimator takes into account all the words of a document by utilizing b n , as well as the length of the document, L n . CHAPTER 3. DATA MINING AND TENSORS

Algorithm 1

 1 Standard averaging for estimating the second order moments Input: D, a corpus of N c documents of length [L 1 , . . . , L Nc ] Output: P Stand-Ave 1: P Stand-Ave = zeros(D,D) 2: for n = 1, 2, . . . , N c do 3:

3 +

 3 ing P and T in (3.15) and (3.16) for Ruffini's estimator are (1+N c ) D(D + 1) 2 + N c and (1 + 2N c ) D+3-1 2N c , respectively. However, these values in standard averaging are N c D(D + 1) 2 + 2N c and N c D+3-1 3 + 2N c . Therefore, in estimating T , standard averaging costs fewer multiplications, and in estimating P , the comparison depends on the values of N c and D. The number of additions in computing P in (3.15) for Ruffini's estimator is (1 + Nc n=1 L n )D + 3N c , while the same quantity for standard averaging is CHAPTER 3. DATA MINING AND TENSORS

Figure 3 . 2 :

 32 Figure 3.2: The consistency of the moment estimates. The relative error in the estimation of second order moments (P ) and third order moments (T ) are plotted in solid and dotted line, respectively. Note that the plot of standard averaging and Ruffini's estimator are superimposed. In this figure, the larger the size of the corpus, the fewer the number of corpuses, to keep the total block length (and hence the computational load) constant.

CHAPTER 3 .

 3 DATA MINING AND TENSORSsider a corpus of N c = 100 documents whose length L n are chosen randomly between a minimum (L min = 3) and a maximum length (L max = 100). The topics and the words of each document are generated according to the generative process described in Section 3.4.1.2. Each experiment averaged over the results of 50 different probabilities φ and A. In these series of experiments, we either plot the moment consistency versus a range of values of D

Figures 3 .

 3 Figures 3.3 and 3.4 show the relative error in estimating P and T versus a range of values of D, when K is fixed to 10. Since in Fig. 3.3, the difference of standard averaging and Ruffini's estimator is not significantly obvious, we plot the results of standard averaging and Ruffini's estimator separately in Fig. 3.4. For the same reason, in the rest of figures of this section, we just show the results of standard averaging and Ruffini's estimator, since the results of simple averaging is much worse than standard averaging and Ruffini's estimator.

Figure 3 . 3 :Figure 3

 333 Figure 3.3: Comparison between all mentioned estimators in estimating P and T with K = 10, D = [10, 20, 30, . . . , 100], N = 100, L min = 3, L max = 100, and averaging over 50 different probabilities φ, A, Left: The relative error of the second order moments estimation, Right: The relative error of the third order moments estimation.

Figure 3 . 5 :

 35 Figure 3.5: Comparison between standard averaging and the Ruffini's estimator in estimating P and T with D = 15, K = [10, 20, 30, . . . , 100], N c = 100, L min = 3, L max = 100, and averaging over 50 different probabilities φ, A, Left: The relative error of the second order moments estimation, Right: The relative error of the third order moments estimation.

Figure 3

 3 Figure 3.6: Comparison between standard averaging and Ruffini's estimator in estimating P and T with D = 65, K = [10, 20, 30, . . . , 100], N c = 100, L min = 3, L max = 100, and averaging over 50 different probabilities φ, A, Left: The relative error of the second order moments estimation, Right: The relative error of the third order moments estimation.

  CHAPTER 4. PERFORMANCE INDEX Section 4.4, and its comparison with other existing measures along with discussions about its advantages are provided in Sections 4.5. The content of this chapter is mainly based on our paper [SCJBZ22], which is published in the journal Signal Processing, Elsevier.

  be the original and estimated matrices respectively, where C M ×N stands for the set of M by N complex-valued matrices. Let us denote the set of permutations of N elements by Perm(N ), and denote by P σ the matrix associated with the permutation σ ∈ Perm(N ). If the columns of A and A are normalized by their L 2 norms, scaling ambiguity reduces to post-multiplication by a diagonal matrix Λ with entries of unit modulus. Assume A = AP σ Λ+W , where the columns of A and A are normalized,

  and denote by C the matrix whose entries are C ij . Then, if the columns of A and A are normalized by their L 2 norms, we have 0 ≤ C ij ≤ 1. In the sequel, three main approaches of measuring the distance between A and A appeared in the literature are reviewed.

  4.3.1 Methods based on correlation matrix4.3.1.1 Greedy approach of [FIW + 20]

  .9) In order to satisfy (4.9) and to consider the scaling ambiguity, a new matrix C with entries C ij = (|a i | -| a j |) 2 is used to define CorrIndex for row vectors

  ∀i and max k C kj = 1, ∀j. Thus CorrIndex(A, A) = 0. Secondly, we prove the converse. If CorrIndex(A, A) = 0, then it implies that max k C ik = 1, ∀i and max k C kj = 1, ∀j. From these two equalities, it can be inferred that there is at least one 1 in each column and row of C. Let us assume C ij = |a H i âj | = 1. According to the Cauchy-Schwarz inequality and the assumption of normalized columns of A and A, we have |a H i âj | ≤ ∥a i ∥∥â j ∥, where the equality of two sides occurs if and only if a i = âj . Since such a conclusion holds for all other associated pairs of columns of A and A, therefore, A = AP σ Λ.

Figure 4

 4 Figure 4.1: CorrIndex and noise. CorrIndex of a random matrix A 6×4 and its permuted noisy version A. This figure confirms the fact that the larger ϵ 0 , the larger CorrIndex.

Figure 4 .

 4 Figure 4.1 confirms the fact that the larger ϵ 0 , the larger CorrIndex.Therefore, in evaluating different decomposition methods, the one with the least CorrIndex would perform the best.In Table4.2, we show the effect of mutual coherence and noise on the results of each index. However, Fig.4.2 investigates the effect of noise intensity, which is measured by SNR. In this experiment, as in Table4.2, we generate a random matrix A of dimension 150×100 with the mutual coherence constant γ = 0.75 and averaged the results over 50 realizations. Then, the matrix A is obtained by permuting randomly the columns of A and by adding a noise matrix according to each SNR value. The goal of this experiment is to show the drawbacks of greedy methods, and to do this, we compare the result of the greedy methods of [FIW + 20, CLDA09, CKAC14] by one of the graphbased method (i.e. MWM[START_REF] Duan | Linear-time approximation for maximum weight matching[END_REF]), which outputs the exact error. Therefore, we can simply conclude the inaccuracy of [FIW + 20, CLDA09, CKAC14].

Figure 4 .

 4 Figure 4.2 shows the relative error (4.13) between A as an estimation of A. As the error output by MWM is exact, the difference between errors output by the greedy methods of [FIW + 20, CLDA09, CKAC14] and the one by MWM show the inaccuracy of [FIW + 20, CLDA09, CKAC14]. As it is expected, the relative error by greedy method [FIW + 20] (resp. [CLDA09, CKAC14]) is optimistic (resp. pessimistic), since its reported error is smaller

Figure 4 . 2 :

 42 Figure 4.2: Drawbacks of greedy methods of [FIW + 20, CLDA09, CKAC14]. Compare the relative error between a random matrix A 150×100 with the mutual coherence constant γ = 0.75 and its permuted noisy version A versus SNR reported by greedy methods of [FIW + 20, CLDA09, CKAC14] and by one of the exact indices, i.e. MWM averaged over 50 realizations.

  .1) holds for every x ∈ R n , where B(x, ϵ) is the closed ball: B(x, ϵ) def = {x|d e (x, x) ≤ ϵ}, (5.2) in which d e (x, x) is the Euclidean distance. See Fig. 5.1 for an example of a lower semi-continuous function. * Sub-gradient [RW09, Definition 8.3]: Suppose R = [-∞, +∞].

f (x) x 0 Figure 5 . 1 :

 051 Figure 5.1: An example of a lower semi-continuous (lsc) function. The solid blue dot indicates f (x 0 ).

CHAPTER 5 .yFigure 5 . 2 :

 552 Figure 5.2: The projection of a vector x ∈ R N onto a closed convex set S ⊂ R N

  5.4)Let Γ 0 (R N ) be the class of lower semi-continuous functions f : R N → (-∞, +∞], with dom(f ) ̸ = ∅. Then i S belongs to Γ 0 (R N ).According to the proposition of Moreau in 1962[START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF], the definition of Proximity operator is obtained by replacing i S (y) in (5.3) with any arbitrary function in Γ 0 (R N ):CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONSProximity operator[START_REF] Patrick | Proximal splitting methods in signal processing[END_REF] For every x ∈ R N , the unique solution of the following minimization problem:

CHAPTER 5 .

 5 CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS can be solved by means of a particular algorithm called Forward-Backward Splitting. Theorem 1 [CP11] explains the relation of minimizing such a cost function with this algorithm.

  as a relaxation parameter, α k , which cannot exceed 1. Algorithm 2 directly DECOMPOSITIONS Algorithm 2 Forward-Backward Splitting [CW05, CP11, Algorithm 10.5] Input: The function f = h + g as defined in Theorem 1, β, x 0 ∈ R N Output: The minimizer of f 1: Fix ϵ ∈ (0, min{1, 1 β })

CHAPTER 5 .

 5 CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS Algorithm 3 Beck-Teboulle proximal gradient algorithm [CP11, Algorithm 10.7] based on FISTA [BT09b] Input: The function f = h + g as defined in Theorem 1, β, x 0 ∈ R N Output: The minimizer of f 1: Set z 0 = x 0 and t 0 = 1 2: for k = 0, 1, 2, . . . do 3:

  g.[START_REF] Comon | Tensors: a brief introduction[END_REF][START_REF] Comon | Independent component analysis, a new concept? Signal Proc[END_REF] for an introduction. Many algorithms have been devised for this kind of decomposition, including the pair-sweeping CoM algorithm, or Joint Approximate Diagonalization (JAD) algorithms[START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF].But authors in [AGH + 14] utilized the tensor power iteration [DLC + 95] to extract the dominant "eigenvector 6 ", â, and the dominant "eigenvalue", φ, of T and then proceeded by deflation, i.e. T ← T -φ â ⊗ â ⊗ â, to get the remaining ones.

  With this in mind, the line 8 of Algorithm 5 includes applying hard-thresholding on each element of Y ([f h √ 2γ (Y (i,j) )] (In,R) (i=1,j=1) ).

5-ℓ 1

 1 , ℓ 2 , ℓ ∞ norms are semi-algebraic. 6-The sum of ℓ 0 pseudo-norm and a polynomial is semi-algebraic [ABS13, Example 5.4]. DECOMPOSITIONS The assumptions in [BST14] are mentioned as several parts in three categories, which are listed below: Assumption 1 (i): r n (.) and r λ (λ) are proper and lsc: since the domains of r n (.) and r λ (λ) are not null, these functions are proper, in addition, as mentioned before, indicator functions and ℓ p norm are lsc according to the definition of lsc on page 70. (ii): f d is a C 1 function (the class of functions with first order differentiability): since f d is quadratic, therefore it is a C 2 function, and hence it is a C 1 function.

9

  We would like to thank the corresponding authors of FastNTF-APG [ZZZ + 16] and BC-VMFB[START_REF] Vu | A proximal approach for nonnegative tensor decomposition[END_REF], Guoxu Zhou and Caroline Chaux, respectively, who sent us the MATLAB codes of their methods. The MATLAB codes of AO-ADMM[START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] and APG[START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF] are made available by the authors at[START_REF] Huang | AO-ADMM-Code @ONLINE[END_REF] and [XY13a], respectively.

Figure 5 . 3 :

 53 Figure 5.3: The reconstruction relative error in decomposing a tensor of dimension 10 × 10 × 10, of rank R = 6 under the non-negativity constraint over all loading factors, in noiseless case with the following setting: ϵ 1 = 10 -20 , average number= 10, initialization number= 10, iterations max-number= 5000, e = 1.9.

Figure 5 . 4 :Figure 5 . 5 :Figure 5 . 6 :Figure 5 . 7 :Figure 5 . 8 :

 5455565758 Figure 5.4: Compare the estimation of matrix X via CorrIndex (Left) and Hungarian (Right) in the same experiment as Fig. 5.3.

Figure 5 . 9 :

 59 Figure 5.9: Compare the estimation of matrix X via CorrIndex (Left) and Hungarian (Right) in the same experiment as Fig. 5.8.

Figure 5 .

 5 Figure 5.12 indicates the relative reconstruction error of a noisy tensor of size 100 × 100 × 100, of rank R = 3 with SNR, initialization number, average number and iterations max-number, equal to 10, 20, 200 and 1000, respectively. As the decomposition of large noisy tensors is a difficult problem to

Figure 5 .

 5 Figure 5.10: The reconstruction relative error in decomposing a tensor of dimension 100 × 100 × 100, of rank R = 3 under the non-negativity constraint over all loading factors, in a noiseless case with the following setting: ϵ 1 = 10 -20 , average number = 10, initialization number = 10, iterations max-number = 5000, e = 1.9.

Figure 5 .Figure 5 . 12 :Figure 5 .Figure 5 .

 551255 Figure 5.11: Compare the estimation of matrix X via CorrIndex (Left) and Hungarian (Right) in the same experiment as Fig. 5.10.

Figure 5 .Figure 5 . 16 :Figure 5 .

 55165 Figure 5.15: Compare the estimation of matrix X via CorrIndex (Left) and Hungarian (Right) in the same experiment as Fig. 5.14.

Figure 5 .Figure 5 .

 55 Figure 5.18: The reconstruction relative error in decomposing a tensor of dimension 10 × 10 × 10, of rank R = 3 under the simplex set constraint over all the columns of all loading factors and over λ, SNR = 10 with the following setting: ϵ 1 = 10 -20 , average number= 200, initialization number= 20, iterations max-number= 1000, e = 1.9.

Figure. 5 .

 5 Figure. 5.20 shows the gap between λ and λ via CorrIndex and Hungarian algorithm (which is the exact error) for the same experiment of Fig. 5.18.

Figure 5 . 19 :Figure 5 .

 5195 Figure 5.19: Compare the estimation of matrix X via CorrIndex (Left) and Hungarian (Right) in the same experiment as Fig. 5.18.

and 5 .

 5 22 correspond to a simulation, which reveals that the estimated probabilities by Robust tensor CHAPTER 5. CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONS

Fig. 5 .

 5 Fig.5.22, we generated some arbitrary probabilities, i.e. A 6×4 and φ 4×1 , which are depicted in Fig.5.21 (left column) and Fig.5.22 (top). Note that in Fig.5.21 (left column), four plots correspond to each column (each hidden variable or topic, i.e. k) of the matrix A, and in each plot, the six values are the conditional probabilities of multi-view variables (words), i.e. d, according to their hidden variable (topic). Similarly, Fig.5.22 (top), contains the probabilities of four hidden variables (topics), i.e. k, in this simulation. Then, by considering these probabilities, synthetic data with hidden relation are generated via the generative procedure described in Section 3.4.1.1 13 . The method of moment estimation in this simulation to obtain empirical second and third order moments is simple averaging with N c = 10 5 documents for averaging. Note that these documents construct a corpus of size N c . Finally, we applied Robust tensor power method and the non-negative AO-ADMM on these empirical moments, and obtained the estimation of A and φ.

Figure 5 . 21 :

 521 Figure 5.21: Left column: an arbitrary conditional probability, A 6×4 . Four plots correspond to each column (each hidden variable or topic, i.e. k) of the matrix A, which contains 6 probability values for each multi-view variables (d). Middle column: the estimation of the matrix A with the non-negative AO-ADMM. Right column: the estimation of the matrix A with Robust tensor power method including some negative estimated probabilities, which are not acceptable.

Figure 5 .

 5 Figure 5.22: Top: an arbitrary probability distribution φ 4×1 contains four values for each hidden variable (topic), k. Middle: the estimation of φ with non-negative AO-ADMM. Bottom: the estimation of φ with Robust tensor power method. Obviously, the result of non-negative AO-ADMM is much more closer than Robust tensor power method to the original probability distribution (top).

Figure 5 .

 5 23 (Left) shows the relative reconstruction error of algorithms versus a range of SNR values, i.e.[10,20, 30, 40]. As it can be seen, Power method is highly sensitive to the additive noise, and its relative reconstruction error can reach very large values (such as 500). In Fig.5.23 (Right), so as to distinguish the performances of other algorithms, we remove Power method.Figure 5.23 (Right)

Figure 5 .Figure 5 .Fig. 5 .

 555 Figure 5.23: Left: The relative reconstruction error of decomposing third order moments including Power method. Right: zoom of the left one, excluding Power method. Averaging the results of decomposing noisy T according to a range of SNR values, i.e. [10, 20, 30, 40] over 200 realizations of A and φ, reveals that constrained algorithms such as SFBS perform much more better than Power method and its variants.

123CHAPTER 5 . 5 (Figure 5 .

 555 Figure 5.25: Moment consistency (cf. Section 3.4.2 for the definition of consistency).

Figure 5 .

 5 Figure 5.26: The comparison of performances of Robust tensor power method and non-negative AO-ADMM in estimating probabilities, i.e. A (top) and φ (bottom), versus corpus size up to N c = 2 17 , when the moment estimator is simple averaging. The number of words and topics are fixed to D = 8 and K = 4, respectively.Each black circle (resp. red cross) corresponds to the error of Robust tensor power method (resp. non-negative AO-ADMM) in estimating φ or A for a specific corpus of documents with a particular corpus size, N c . The median (resp. standard deviation) among corpuses of documents is also plotted in solid (resp. dotted) line for every corpus size.

Moreover, Fig. 5 .

 5 26 (bottom) reveals that the difference in the performance and standard deviation of Robust tensor power method is much larger in estimating φ. Therefore, taking into account a proper constraint (such as non-negativity) in the decomposition is very effective, especially in the estimation of φ.Results for random A and φ. When running extensive computer simulations, one can assume a uniform distribution for the hidden variables (topics), and K uniform distributions for the conditional probability of multiview variables (words) given the topic, which are actually the columns of A. Once A and φ are obtained, one can proceed exactly the same way as explained in Subsection 3.4.1 to generate synthetic datasets and compute sample moments via simple averaging.

Figure 5 .

 5 Figure 5.27 (top) compares the performances of algorithms in estimatingA versus the corpus size. As it can be seen in Fig.5.27 (top), the performances of all algorithms except Power method are the same, which shows that the projection to the simplex set or considering this constraint in the

Figure 5 .

 5 Figure 5.27: The comparison of performances in estimating probabilities, i.e. A (top) and φ (bottom), versus corpus size up to N c = 2 10 , when the moment estimator is standard averaging. The number of multi-view variables (words) and the number of hidden variables (topics) are fixed to D = 10 and K = 3, respectively.Robust tensor power method is run with M 1 = 5000 and M 2 = 20, and the initialization number and iterations max-number for constrained algorithms is 20 and 1000, respectively.

CHAPTER 5 .

 5 CONSTRAINED AND NON-CONSTRAINED TENSOR DECOMPOSITIONSdecomposition will end in more precise estimations compared to pure Power method.

Figure 5 .

 5 Figure 5.27 (bottom) compares the performances of algorithms in estimating φ versus the corpus size. Unlike the estimation of A, the performance of different algorithms are not the same. Actually, the performance of Projected Power method is the worst, while the relative error of φ by Symmetric Simplex SFBS is the least among other methods. Therefore, Symmetric Simplex SFBS outperforms other methods in decomposing the tensor of third order moments.

B. 1 Figure B. 1 :Algorithm 6

 116 Figure B.1: Le modèle à sujet unique avec ses variables multi-vues (x i ) et leur variable cachée correspondante (h).

12 :

 12 for toutes les paires de mots possibles, ex. (word 1 , word 2 ) do 7:P temp (word 1 , word 2 ) = 1 + P temp (word 1 , word 2 ) 8:count-pair = 1 + count-pair Symétriser P Stand-Ave de la taille du dictionnaire (D), ce qui est un avantage dans le cas de dictionnaires de grande dimension.

B. 3 |S|S|S

 3 .1.2 Méthodes basées sur le transport optimal Problème d'affectation optimale Le problème d'assignation est un problème d'optimisation combinatoire ancien, bien connu et fondamental [Gal83, Mun57, PC + 19]. Le premier algorithme en temps polynomial pour les problèmes d'affectation optimale est la "méthode hongroise" [Kuh55] également connue sous le nom de "Kuhn-Munkres" [Mun57,TK04], et la complexité de l'algorithme est approximativement O(N 4 ) flops [Mun57]. Cet algorithme a été utilisé dans [TK04] pour un appariement optimal des sources dans BSS. B.3.1.3 Indices invariants à la permutation Cependant, dans la littérature sur la séparation des sources [CJ10], certains indices ont été proposés pour mesurer l'écart (basé sur une définition spécifique de l'écart) entre les matrices de mélange originales et estimées sans chercher à trouver la permutation correspondante [Com94,MM94,ACY + 96]. Les indices proposés dans [Com94, MM94, ACY + 96] sont basés sur S = A -1 A (ou S = A † A pour les matrices A non carrées). Les détails de ces indices sont les suivants. Indice de Comon [Com94] L'indice de Comon est une combinaison des normes L 1 et L 2 , et se calcule comme suit : ij | 2 -1 . La valeur de ϵ 1 peut augmenter énormément, en fonction des valeurs de la matrice S, par conséquent, cet indice n'est pas borné supérieurement. Indice de Moreau-Macchi [MM94] L'indice proposé dans [MM94] mesure un écart entre la matrice S et une matrice de permutation. Il est défini comme suit :ij | 2 (max k |S ik |) 2 -ij | 2 (max k |S kj |) 2 -1 .La division par la valeur maximale (e.g. (max k |S ik |) 2 ) fournit une borne supérieure pour ϵ 2 contrairement à ϵ 1 . Indice de Amari [ACY + 96] Cet indice de performance prend la forme : max k |S kj | -1 .

Figure B. 5 :

 5 Figure B.5: CorrIndex et bruit. CorrIndex d'une matrice aléatoire A 6×4 et de sa version bruitée permutée A. Cette figure confirme le fait que plus le bruit est grand, plus le CorrIndex est grand.

Theorem 3 (

 3 Séparation Avant-Arrière[START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]). Supposons que f : R N → R ∪ {+∞} est une fonction propre 3 , semi-continue inférieurement, qui possède la propriété KL et est bornée par le bas. Si f peut être décomposée en deux parties sous la forme f = h + g, où g est semi-continue inférieurement et h : R N → R est une fonction différentiable à valeur finie avec un gradient continu β-Lipschitz, i.e, ∃β tel que :∥∇h(x) -∇h(y)∥ 2 ≤ β∥x -y∥ 2 ,alors on peut montrer[START_REF] Patrick | Signal recovery by proximal forward-backward splitting[END_REF] que le minimiseur de f satisfait l'équation de point fixe suivante :x = prox γg (x -γ∇h(x)), (B.19)où γ ∈ (0, +∞).Plusieurs variantes d'implantation de l'agoritme de Séparation Avant-Arrière existent, et sont rapportées dans[START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]. Deux d'entre elles sont reprises (Algorithme 7 et Algorithme 8) auxquelles nous ferons référence dans la suite de cette thèse.L'algorithme 8 est basé sur l'algorithme Fast Iterative Shrinkage Thresholding (FISTA) proposé dans[START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF][START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF] et peut être considéré comme un algorithme de gradient proximal. En fait, au lieu de la fonction de rétrécissement dans FISTA, on utilise l'opérateur de proximité.Algorithm 7 Séparation Avant-Arrière [CW05, CP11, Algorithm 10.5] Entrée: La fonction f = h + g telle que définie dans le Théorème 1, β,x 0 ∈ R N Sortie: Le minimiseur de f 1: Fixer ϵ ∈ (0k = x k -γ k ∇h(x k ) 5: α k ∈ [ϵ, 1]6:x k+1 = x k + α k (prox γ k g (y k ) -x k )7: end for La convergence de l'algorithme de Séparation Avant-Arrière dans le théorème 1 est prouvée dans [ABS13], et l'analyse de convergence fournie n'est pas seulement applicable pour des fonctions de coût continues et convexes, mais est également utilisable pour des fonctions non-lisses et non-convexes (ou un ensemble de contraintes non-convexes). B.4.3 L'état de l'art B.4.3.1 Optimisation Alternée -Méthode des Multiplicateurs à Direction Alternée (AO-ADMM) [HSL16] L'algorithme AO-ADMM 4 tente de minimiser chaque étape de (5.9) par ADMM. Les détails de l'AO-ADMM pour la décomposition CP sous contrainte peuvent être trouvés dans [HSL16]. La convergence d'AO est brièvement examinée dans [HSL16] sur la base des travaux de [Tse01, RHL13]. La convergence de chaque étape, qui est un algorithme ADMM, peut être démontrée pour les fonctions convexes [HSL16, BPC + 11] et aussi récemment pour certaines fonctions non-convexes telles que ℓ q , 0 < q < 1 [WYZ19], mais pas encore pour certaines autres telles que ℓ 0 . 4 Alternating Optimization -Alternating Direction Method of Multipliers 168 APPENDIX B. RÉSUMÉ EN FRANCAIS Algorithm 8 Algorithme du gradient proximal de Beck-Teboulle [CP11, Algorithm 10.7] basé sur FISTA [BT09b] Entrée: La fonction f = h + g telle que définie dans le Théorème 1, β, x 0 ∈ R N Output: Le minimiseur de f 1: Set z 0 = x 0 and t 0 = 1 2: for k = 0, 1, 2, . . . do 3:

) 8 : 5 :

 85 end for B.4.3.2 Gradient Proximal Alternatif (APG)[START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF] Dans[START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF], la convergence d'un algorithme BCD général est étudiée, où pour mettre à jour chaque bloc de variables inconnues, trois types de mises à jour sont étudiés, à savoir : original, proximal et prox-linear. D'après la comparaison effectuée dans [XY13b], la mise à jour proxlinéaire (i) donne de meilleures valeurs de la fonction objectif que la mise à jour originale et proximale, (ii) est plus facile à calculer et (iii) permet souvent des solutions sous forme explicite. L'algorithme APG 5 ne suit pas explicitement la procédure de Séparation Avant-Arrière; néanmoins, pour satisfaire la contrainte de non-négativité, il se termine par une projection sur l'orthant non-négatif ainsi que par la mise à jour d'un coefficient particulier comme à la ligne 5 de l'Algorithme 8. Algorithm 9 L'algorithme de SFBS Entrée: T , C A (n) , initiale A (n) 0 , n ∈ [1, . . . , N ], e Sortie: Estimation de A (n) , n ∈ [1, . . . , N ] 1: repeat 2: for n = 1, 2, . . . , N do 3: W = (A (N ) ⊙ . . . ⊙ A (n+1) ⊙ A (n-1) ⊙ . . . ⊙ A (1) ) T 4: β = {max(valeur singulière(W ))} 2 définir γ = e β et choisir α g .6:

LaFigure B. 7 :

 7 Figure B.7: Erreur relative de reconstruction lors de la décomposition d'un tenseur de dimensions 10 × 10 × 10, de rang R = 6 sous la contrainte de non-négativité sur tous les matrices facteurs, dans le cas sans bruit avec les paramètres suivants: ϵ 1 = 10 -20 , nombre moyen= 10, nombre d'initialisation= 10, itérations max-number= 5000, e = 1.9.

Figure B. 8 :Figure B. 9 :

 89 Figure B.8: Comparez l'estimation de la matrice X via CorrIndex (à gauche) et Hungarian (à droite) dans la même expérience que Fig. B.7.

Figure B. 10 :Figure B. 11 :

 1011 Figure B.10: Comparaison entre l'estimation de la matrice X via CorrIndex (à gauche) et Hungarian (à droite) dans la même expérience que la Fig. B.9.

Figure B. 12 :

 12 Figure B.12: Comparaison entre l'estimation de la matrice X via CorrIndex (à gauche) et Hungarian (à droite) dans la même expérience que la Fig.B.11.

Figure B. 14 :Figure B. 15 :

 1415 Figure B.14: Comparaisonentre l'estimation de la matrice X via CorrIndex (à gauche) et Hungarian (à droite) dans la même expérience que la Fig. B.13.

  Notation Throughout the thesis, arrays of numbers will be printed in boldface. More precisely, one-way and two-way arrays will be denoted in bold lowercase and bold uppercase, respectively, e.g. v and M . Arrays with more than two indices will be denoted by bold calligraphic symbols, such as

	CHAPTER 1. INTRODUCTION
	are devoted to the main contributions of this thesis, describe performance in-
	dices (including CorrIndex) and tensor decomposition algorithms (including
	SFBS), respectively. Finally, Chapter 6 concludes the thesis and mentions
	some perspectives.
	2 Tensor

Chapter 3, explains the relation between tensor decompositions and the data mining problem, which is in the focus of attention of this thesis. Chapter 4 and Chapter 5, which A. Sets and spaces will be noted in script font, like S. Entries of arrays v, M and A will be noted v i , M ij and A ijk , without bold font (since they are indeed scalar numbers).

Contents 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . 2.2 Tensor: notations and preliminaries . . . . . . . 2.3 Exact tensor decomposition . . . . . . . . . . . . 2.3.1 Canonical Polyadic (CP) tensor decomposition . . 2.3.2 Tucker tensor decomposition . . . . . . . . . . . . 2.3.3 Multi-Linear/Higher Order Singular Value Decomposition (MLSVD/HOSVD) . . . . . . . . . . . . . 2.4 Approximate tensor decomposition . . . . . . . . 2.4.1 Truncated MLSVD/HOSVD . . . . . . . . . . . . 2.4.2 Low-rank CP approximation . . . . . . . . . . . . 2.5 Constrained CP decomposition . . . . . . . . . . 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . .

Table 3

 3

.2: Compare the number of required multiplications Estimator

Table 4 .

 4 1: Approximate numbers of multiplications of computing each stage of

	CorrIndex and other methods
	Method	C or S Permutation Index	Significance
			of upper bound

Table 4 .

 4 2: A numerical comparison on methods measuring the distance betweenA 150×100 with the mutual coherence constant γ = 0.75 and its permuted noisy version A with SNR = -1.76 dB averaged over 50 realizations. The index in five first rows of the table is the relative error defined in (4.13). On the other hand, the last four indices of the table are defined differently and are hence not comparable.

	Method	Index	Computing	Significance
			time (ms)	of upper bound
	Greedy of [FIW + 20]	0.37	0.9	No
	Greedy of [CLDA09, CKAC14] 1.05	5.3	No
	Hungarian [Mun57]	0.86	4.5	Yes
	MWM [DP14]	0.86	2.9	Yes
	Linprog [PC + 19]	0.86	1310	Yes
	Comon [Com94]	1.8e4	3.3	No
	Moreau-Macchi [MM94]	897.91	3.3	No
	Amari [ACY + 96]	3.2e3	2.9	No
	CorrIndex	0.36	0.4	Yes

Table 4 .

 4 3: A numerical comparison on methods measuring the distance betweenA 150×100 with the mutual coherence constant γ = 0.95 and its permuted noiseless version A averaged over 50 realizations. The index in five first rows of the table is the relative error defined in (4.13). The four last indices of the table are defined differently and are hence not comparable.The experiment of Table 4.2 corresponds to an inaccurate estimation ofA (i.e. SNR = -1.76), and all the performance indices (perhaps except the greedy method of [FIW + 20], which yields 37%) demonstrate convincingly the fact that estimation is not accurate. In order to evaluate the indices in the opposite situation (accurate estimation), we perform another experiment with the same setting as that of Table 4.2 except that A is a permuted noiseless version of A 150×100 with γ = 0.95. The result of this experiment is reported in Table 4.3.

	Method	Index	Computing	Significance
			time (ms)	of upper bound
	Greedy of [FIW + 20]	3.82e -17	0.89	No
	Greedy of [CLDA09, CKAC14]	0.10	5.2	No
	Hungarian [Mun57]	0	2.7	Yes
	MWM [DP14]	0	4.0	Yes
	Linprog [PC + 19]	0	1630	Yes
	Comon [Com94]	2.14e -13	3.3	No
	Moreau-Macchi [MM94]	0	4.5	No
	Amari [ACY + 96]	3.35e -16	3.8	No
	CorrIndex	2.37e -16	0.59	Yes
	different from (4.13).

As can be seen in Table

4

.2, CorrIndex is the fastest index. In addition, contrary to the indices of Comon, Moreau-Macchi and Amari, Cor-rIndex returns a value in the bounded range [0, 1]. One could normalize the Moreau-Macchi and Amari indices in order to obtain bounded values, but the signification of such upper bounds on these indices has never been investigated.

Table 5 .

 5 1: Proximity operator of functions used in this thesis

	Function	Proximity operator
	Indicator function of S (i S (5.4)) Projection onto S (proj S (5.3))
	ℓ 0 pseudo-norm (cardinality)	hard-thresholding (5.14)
	ℓ 1 norm	soft-thresholding (5.16)
	5.5.2.4 Sparsity with ℓ 1 norm

A common and convex approximation of the ℓ 0 pseudo-norm is ℓ 1 norm.

The resulting cost function obtained by replacing ℓ 0 in (5.15) with ℓ 1 norm is called LASSO regression

[START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] David | Compressed sensing[END_REF]

.

Table 5 .

 5 2: Generating some arbitrary matrices A 10×3 and vectors φ 3×1 , then calculating their corresponding T and P , averaging the results of decomposing noiseless T over 200 realizations of A and φ.

	Algorithm	Reconstruction Error of A Error of φ
	Power method	2.97 e-15	5.25 e-16 3.68 e-16
	Projected Power method	9.19 e-16	5.10 e-16 7.70 e-16
	Simplex AO-ADMM	4.11 e-5	1.50 e-4	2.40 e-4
	Simplex SFBS	1.57 e-15	4.97 e-15 8.39 e-15
	Symmetric Simplex SFBS	1.52 e-15	4.97 e-15 8.39 e-15

Table 5 .

 5 3: The median and standard deviation of error in estimating A and φ, with a corpus of size N c = 2 17 ≈ 1.2 × 10 6 , and K = 4. In each cell, top: non-negative AO-ADMM, bottom: Robust tensor power method. As it can be seen, the median and standard deviation of the error of non-negative AO-ADMM is always less than those of Robust tensor power method.

	Median of relative error in estimating A
	D	6	7	8	9	10
	Non-negative AO-ADMM	0.09 0.08 0.07 0.10 0.07
	Robust tensor power method 0.38 0.18 0.24 0.18 0.12
	Standard deviation of relative error in estimating A
	D	6	7	8	9	10
	Non-negative AO-ADMM	0.10 0.12 0.08 0.10 0.12
	Robust tensor power method 0.31 0.31 0.39 0.33 0.23
	Median of relative error in estimating φ
	D	6	7	8	9	10
	Non-negative AO-ADMM	0.14 0.10 0.11 0.10 0.13
	Robust tensor power method 0.66 0.28 0.31 0.29 0.23
	Standard deviation of relative error in estimating φ
	D	6	7	8	9	10
	Non-negative AO-ADMM	0.22 0.12 0.21 0.28 0.31
	Robust tensor power method 0.54 0.52 0.59 0.41 0.40

Table 5 .

 5 4: Performances in estimating probabilities φ and A in terms of CorrIndex:

	the smaller the better. This experiment is carried out using a part of a well-known
	text data set, namely 20 Newsgroups. To be more precise, documents of the selected
	data set belong to four topics: "computer graphics", "baseball", "cryptography" and
	"Christianity".
	Algorithm	CorrIndex(φ, φ) CorrIndex(A, A)
	Power method	0.577	0.077
	Projected Power method	0.703	0.077
	Simplex SFBS	0.106	0.097
	Symmetric Simplex SFBS	0.106	0.097
	Simplex AO-ADMM	0.102	0.089

  Table B.1: Comparaison numérique des méthodes de mesure de la distance entre A 150×100 avec une cohérence mutuelle γ = 0, 75 et sa version bruitée permutée A avec SNR = -1, 76 dB en moyenne sur 50 de réalisations. L'indice dans les cinq premières lignes du tableau est l'erreur relative. En revanche, les quatre derniers indices du tableau sont définis différemment et ne sont donc pas comparables.

	Méthode	Index	Temps de	Significativité
			calcul (ms)	de la limite supérieure
	Greedy of [FIW + 20]	0.37	0.9	Non
	Greedy of [CLDA09, CKAC14] 1.05	5.3	Non
	Hungarian [Mun57]	0.86	4.5	Oui
	MWM [DP14]	0.86	2.9	Oui
	Linprog [PC + 19]	0.86	1310	Oui
	Comon [Com94]	1.8e4	3.3	Non
	Moreau-Macchi [MM94]	897.91	3.3	Non
	Amari [ACY + 96]	3.2e3	2.9	Non
	CorrIndex	0.36	0.4	Oui

qui donnent lieu à des estimations plus fiables et plus précises. De plus, nous introduisons un algorithme de décomposition tensorielle sous contraintes, appelé Simple Forward-Backward Splitting (SFBS), qui est basé sur une approche de minimisation proximale. SFBS est plus performant que l'état de l'art dans la décomposition des tenseurs bruités tout en étant moins coûteux en calcul.En outre, pour évaluer la performance des algorithmes de décomposition tensorielle, nous introduisons un indice appelé CorrIndex, qui fournit des limites de performance interprétables, tout en maintenant la complexité de calcul à un niveau raisonnable. De plus, nous proposons une méthode d'estimation des moments (moyenne standard), qui estime les moments d'ordre deux et trois, avec les mêmes performances que l'état de l'art, mais basée sur un concept beaucoup plus simple, i.e. la moyenne pondérée. De plus, le calcul de la moyenne standard est plus performant dans les petites dimensions et présente certains avantages en termes de complexité de calcul.

The Kruskal rank of a matrix is defined as the maximum number k such that any k columns of that matrix are linearly independent.

This can be done with the commands mnrnd and numpy.random.multinomial in MAT-LAB and Python, respectively.

convex by definition.

For two arbitrary functions p(.) and q(.), we say p(x) = O(q(x)) as x → x, if there exist positive values δ and M such that for all 0 < |x -x| < δ, |p(x)| < M q(x).

A differentiable function has a single sub-gradient at each point.

A function is proper, if its domain is not a null set.

Subsequence converging to a Nash point[START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF] Equation 2.3]

The spectral norm of a matrix is defined as its maximum singular value[START_REF] Bernstein | Matrix Mathematics[END_REF].

Being finite-value is not boundedness, but it means that a function takes values in the real line, i.e. (-∞, +∞).

By hidden relation, we mean the relation between hidden/latent variables (topics) and multi-view variables (words), which is depicted in Fig.3.1.

For example, if we report ϵ( T ) = 1.5, it should not be interpreted as 1.5%, but as 150%.

In plots, we show the loading factors of a third order tensor by A, B, C.

To calculate cumulative distributions, the MATLAB function cumsum() can be used.

In this experiment, we work with true moments instead of their estimations. Therefore, we do not need realizations of documents. Note that we need realizations of A and φ to generate different values of T and P , but with realizations of documents, the empirical estimations in (3.8) and (3.7) are computed.

We did these steps by means of corresponding packages such as "nltk", libraries such as "Gensim" and an implemented example available online[START_REF] Prabhakaran | 20 Newsgroups -Pre-processing steps @ONLINE[END_REF].

Une fonction est propre si son domaine n'est pas un ensemble nul.

Alternating Proximal Gradient

Voir page 101 et Algorithme 5 pour la définition de ces paramètres.

Scientific Computing, 78(1):29-63, 2019.

Acknowledgments

CHAPTER 4. PERFORMANCE INDEX

Recently, by improving some required computational steps, the running time of implementing linear programming algorithms has been reduced to approximately q 2+ 1 6 [CLS21] and q 2+ 1 18 [START_REF] Tat | Solving empirical risk minimization in the current matrix multiplication time[END_REF] flops, where q is the size of unknown vector in the linear programming problem. Therefore, as q = N 2 in (4.7), the lowest complexity to find the optimal permutation of the problem described at the beginning of Section 4.3 by the means of linear programming is approximately N 4 flops. In addition, in [PC + 19, KA21], authors proposed sub-optimal solutions, called Auction algorithms, which costs approximately N 2 flops [START_REF] Khosla | Revisiting the auction algorithm for weighted bipartite perfect matchings[END_REF].

Indices invariant to permutation

Methods described in Section 4.3.1 try first to estimate the permutation, and then measure some distances based on the estimated permutation. As it is seen in Section 4.3.1.1, these methods may actually not return a permutation, and there is no guarantee that ϵ 0 (A, A) = 0 even if the indices they output are zero. Conversely, the algorithms of Section 4.3.2 behave better (in terms of returning optimal permutation) but may become very costly for large values of N .

However, in the literature of source separation [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF], some indices have been proposed to measure the gap (based on a specific definition of gap) between original and estimated mixing matrices without inquiring to find the corresponding permutation [Com94, MM94, ACY + 96]. Moreover, these indices are zero if and only if ϵ 0 (A, A) = 0, which offers a valuable guarantee.

The indices proposed in [Com94, MM94, ACY + 96] are based on S = A -1 A (or S = A † A for non-square A). The details of these indices are as follows.

CHAPTER 4. PERFORMANCE INDEX

for instance in tensor decompositions or in blind source separation. Existing performance indices are classified in three main categories: "greedy methods", "graph-based methods" and "invariant indices". These methods are reviewed, and it is inferred that greedy methods are not reliable especially in noisy situations (they are either optimistic or pessimistic). In addition, graphbased methods and invariant indices are computationally expensive. We propose a new performance index belonging to the class of invariant indices, named CorrIndex, whose upper and lower bounds are easy to interpret, while being computationally cheap. In the rest of this thesis, we will evaluate the performance of tensor decomposition methods by means of CorrIndex. DECOMPOSITIONS Since the relative error of loading factors (like what are represented in Fig. 5.5 -Fig. 5.7) results almost the same result as the gap between X and X via CorrIndex and Hungarian, in the rest of the experiments of this section, we shall just report the relative reconstruction error and the gap between X and X via CorrIndex and Hungarian (like Fig. 5.4). As in Fig. 5.3, BC-VMFB does not perform properly in the experiment of Fig. 5.8. In addition, although Nway carries out well in a noiseless case, its performance is not acceptable in the noisy situation (cf. Fig. 5.8). As mentioned before, FastNTF-APG replaces the noisy tensor with its low rank approximation, which helps to filter the noise out. As it can be seen in Fig. 5.8, the result of FastNTF-APG is better than APG in a noisy scenario, whereas APG is one the best methods for the decomposition of a noiseless tensor (cf. Fig. 5.3). Moreover, SFBS and AO-ADMM outperform other methods in the experiment of Fig. 5.8. Therefore, the only algorithm that performs properly and better than others in both noiseless and noisy situations is SFBS.

In Fig. 5.9, the gap between X and X via CorrIndex and Hungarian algorithms (which is the exact error) is reported for the same experiment of Fig. 5.8. Although according to the relative reconstruction error in Fig. 5.8, both AO-ADMM and SFBS outperform others, in estimating loading factors (cf. Fig. 5.9), SFBS performs better than AO-ADMM. The performance on estimating loading factors is critical in many applications such as data mining and text mining explained in Chapter 3. In addition, as Nway performs inadequately in the noisy case, this reveals the importance of considering constraints in the decomposition of tensors (remind that Nway executes simply non-constrained CP decomposition). DECOMPOSITIONS simplex constraint is directly considered and investigated by the authors in [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] as well as in their implementation. Secondly, according to the results under the non-negativity constraint in Section 5.6.1.1, AO-ADMM is the only algorithm, which has the closest performance to that of SFBS in both noiseless and noisy cases.

The methodology for applying the simplex set constraint on all the columns of all loading factors in addition to the vector of coefficients (λ) is explained in Section 5.5.2.2, and we have utilized this methodology in our simulations of the current section. Noiseless cases. Figure 5.14 shows the relative reconstruction error of the noiseless tensor of size 10 × 10 × 10, of rank R = 3 with initialization number, average number, e and iterations max-number, equal to 10, 10, 1.5 and 20000, respectively, under the simplex set constraint over all the columns of all loading factors in addition to λ 12 . By comparing Fig. 5.14 and Fig. 5.3, it can be inferred that the simplex set is a more difficult constraint than non-negativity to be achieved, since the maximum required iterations for the simplex set constraint is 20000, and it is more than what is set for the non-negativity constraint (i.e. 5000) for the convergence of the algorithms.

As it can be seen in Fig. 5.14, both SFBS and AO-ADMM achieve the same level of relative error, however SFBS converges slightly faster.

In Fig. 5.15, the gap between X and X via CorrIndex and Hungarian algorithm (which is the exact error) is reported for the same experiment of Fig. 5.14. The same conclusion as before can be drawn from estimating loading factors, that is, both SFBS and AO-ADMM perform well, however, SFBS converges a bit faster.

In order to show that Fig. 5.15 has sufficient information even under the simplex set constraint (like the non-negativity constraint), we bring the result of the estimation of first mode loading factor (A) in Fig. 5.16, which has the same information as Fig. 5.15. As matrix X (cf. Fig. 5.15) does not 12 See page 101 and Algorithm 5 for the definition of these parameters.

APPENDIX A. LIPSCHITZ CONSTANT OF THE GRADIENT OF THE FIDELITY TERM IN (??)

Note that the definition of "spectral norm" of a matrix is:

which is equal to the maximum singular value of A with X as a matrix.

Therefore, we have:

Soit L le nombre de mots dans un document donné, x ℓ les mots observés, Comme mentionné précédemment, l'objectif principal est d'estimer la probabilité de la variable cachée (sujet) et la probabilité conditionnelle des variables multi-vues (mots).

À cette fin, certaines hypothèses sont nécessaires et sont énumérées cidessous [AGH + 14, AHK12a] :

• Les probabilités conditionnelles ne dépendent pas de l'ordre des mots (échangeabilité), e.g.

Prob(u γe(p) , u γe(q) , u γe(r) |h) = Prob(u γe(q) , u γe(r) , u γe(p) |h)

• Les mots du sujet sont conditionnellement indépendants compte tenu du sujet, i.e.

(u γe(p) |h) ⊥ ⊥ (u γe(q) |h)

• Les mots du sujet ont la même distribution conditionnelle étant donné le sujet, i.e.

(u γe(p) |h) ∼ (u γe(q) |h).

où ⊗ désigne le produit tensoriel, P est une matrice symétrique D × D et

En raison de ce choix pour u d , ces moments présentent les relations suivantes :

B.3 Indices de performance

Les ambiguïtés de permutation et d'échelle sont des problèmes pertinents dans des applications telles que la décomposition tensorielle [START_REF] Comon | Tensors: a brief introduction[END_REF] et la séparation aveugle de sources (BSS) [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]. La mise à l'échelle fait référence à la multiplication par une matrice diagonale à composantes non nulles, qui peut être complexe dans le cas le plus général, et la permutation fait référence à la permutation des colonnes d'une matrice, ce qui est équivalent à la multiplication par une matrice de permutation. (B.12) L'affectation proposée par cette méthode est (â 1 , a 2 ), (â 2 , a 2 ), (â 3 , a 3 ), ce qui n'est évidemment pas acceptable car la colonne a 2 est sélectionnée deux fois. Si on calcule l'erreur quadratique via 1 2 3 n=1 ∥a n -âσ(n) ∥ 2 2 en considérant l'hypothèse de normalisation de a n et âσ(n) par rapport à la norme L 2 Matching pondéré maximum (MWM 1 ) Le problème de l'affectation optimale peut également être considéré comme un cas particulier de l'matching pondéré maximum (MWM), qui est un problème bien connu en théorie des graphes, pour lequel plusieurs algorithmes en temps polynomial existent [START_REF] Galil | Efficient algorithms for finding maximal matching in graphs[END_REF]. Les meilleurs algorithmes de MWM exacts [DP14,DK69,D + 59] et approximatifs [START_REF] Harold | Faster scaling algorithms for network problems[END_REF] coûtent environ 8N 3 et N 2 flops, respectivement.

Programmation linéaire La recherche de la permutation optimale σ, i.e pour la matrice de permutation optimale P ⋆ , peut être formulée ceci comme le problème d'optimisation suivant [PC + 19] :

Récemment, en améliorant certaines étapes de calcul requises, le temps d'exécution des algorithmes de programmation linéaire a été réduit à environ

flops, où q est la taille du vecteur inconnu dans le problème de programmation linéaire. Par conséquent, comme q = N 2 dans (B.14), la complexité la plus faible pour trouver la permutation optimale du problème décrit au début de la section 4.3 au moyen de la programmation linéaire est approximativement de N 4 flops.

n 3 ×R , sous une forme dépliante, sont les suivants :

k )∥ F , mais comme nous le prouvons dans l'annexe A, la constante de Lipschitz est en fait la norme spectrale 7 de {(A Ensuite, cette matrice de blanchiment est appliquée au tenseur T pour donner :

En d'autres termes, les composantes du tenseur blanchi sont, d'après (2.2) [START_REF] Comon | Tensors: a brief introduction[END_REF]:

Ce nouveau tenseur satisfait

Un problème général de décomposition CP contrainte de T peut être formulé comme suit :

.

Alors en développant en mode-n l'expression (5.11), on a : Bien que dans ce corpus, le nombre de documents pour tous les sujets soit presque le même (ce qui signifie que φ devrait avoir une distribution proche de la distribution uniforme), le vecteur φ estimé par la méthode de puissance φ Power-proj = 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0

B.5 Conclusion et perspectives

Dans cette thèse, nous nous sommes concentrés sur un modèle de données particulier, appelé single topic (cf.