Rima Khouja

Pierre-Alexandre Mattei

Bernard Mourrain

Tensor decomposition for learning Gaussian mixtures from moments

In data processing and machine learning, an important challenge is to recover and exploit models that can represent accurately the data. We consider the problem of recovering Gaussian mixture models from datasets. We investigate symmetric tensor decomposition methods for tackling this problem, where the tensor is built from empirical moments of the data distribution. We consider identifiable tensors, which have a unique decomposition, showing that moment tensors built from spherical Gaussian mixtures have this property. We prove that symmetric tensors with interpolation degree strictly less than half their order are identifiable and we present an algorithm, based on simple linear algebra operations, to compute their decomposition. Illustrative experimentations show the impact of the tensor decomposition method for recovering Gaussian mixtures, in comparison with other state-of-the-art approaches.

Introduction

With the relatively recent evolutions of information systems over the last decades, many observations, measurements, data are nowadays available on a variety of subjects. However, too much information can kill the information and one of the main challenges remains to analyse and to model these data, in order to recover and exploit hidden structures.

To tackle this challenge, popular Machine Learning technologies have been developed and used successfully in several application domains (e.g. in image recognition [START_REF] He | Deep residual learning for image recognition[END_REF]). These techniques can be grouped in two main classes: Supervised machine learning techniques are approximating a model by optimising the parameters of an enough general model (e.g. a Convolution Neural Network) from training data. Unsupervised machine learning techniques are deducing the parameters characterising a model directly from the given data, using an apriori knowledge on the model. The supervised approach requires annotated data, with a training step that can introduce some bias in the learned model. The unsupervised approach can be applied directly on a given data set avoiding the costly step of annotating data, but the quality of the output strongly depends on the type of models to be recovered.

We consider the latter approach and show how methods from effective algebraic geometry help finding hidden structure in data that can be modelled by mixtures of Gaussian distributions. The algebraic-geometric tool that we consider is tensor decomposition. It consists in decomposing a tensor into a minimal sum of rank-1 tensors. This decomposition generalises the rank decomposition of a matrix, with specific and interesting features. Contrarily to matrix rank decomposition, the decomposition of a tensor is usually unique (up to permutations) when the rank of the tensor, that is the minimal number of rank-1 terms in a decomposition, is small compared to the dimension of the space(s) associated to the tensor. Such a tensor is called identifiable. This property is of particular importance when the decomposition is used to recover the parameters of a model. It guaranties the validity of the recovering process and its convergence when the number of data increases.

It has been shown in [START_REF] Chiantini | On generic identifiability of symmetric tensors of subgeneric rank[END_REF] that for symmetric tensors, if the rank of the tensor is strictly less than the rank r g of a generic tensor of the same size, then the tensor is generically identifiable, except in three cases. We show in Theorem 3.6 a more specific result: for a symmetric tensor T having a decomposition with r points, if the Hankel matrix associated to T in a degree strictly bigger than the degree of interpolation of the r points is of rank r, then the tensor is identifiable. We show in Proposition 3.3, that under some assumption on the spherical gaussian mixtures, a tensor of moments of order 3 of the distribution is identifiable and its decomposition allows to recover the parameters of the Gaussian mixture.

Several types of method have been developed to tackle the difficult problem of tensor decomposition. Direct methods based on simultaneous diagonalisation of matrices built from slices of tensors have been investigated for 3rd order multilinear tensors, e.g. in [START_REF] Harshman | Foundations of the parafac procedure: Models and conditions for an "explanatory" multi-modal factor analysis[END_REF][START_REF] Sanchez | Tensorial resolution: A direct trilinear decomposition[END_REF][START_REF] Leurgans | A Decomposition for Three-Way Arrays[END_REF][START_REF] Domanov | Canonical Polyadic Decomposition of Third-Order Tensors: Reduction to Generalized Eigenvalue Decomposition[END_REF] or for multilinear tensors of rank smaller than the lowest dimension in [START_REF] De | A Link between the Canonical Decomposition in Multilinear Algebra and Simultaneous Matrix Diagonalization[END_REF][START_REF] Luciani | Canonical Polyadic Decomposition based on joint eigenvalue decomposition[END_REF]. In his proof on lower bounds of tensor ranks, Strassen showed in [START_REF] Volker Strassen | Rank and optimal computation of generic tensors[END_REF]Theorem 4.1] that a 3rd order multilinear tensor is of rank r if it can be embedded into a tensor with slices of rank r matrices, which are simultaneously diagonalised.

For symmetric tensor decomposition, a method based on flat extension of Hankel matrices or commutation of multiplication operators has been proposed in [START_REF] Brachat | Symmetric tensor decomposition[END_REF] and extended to multi-symmetric tensors in [START_REF] Bernardi | General tensor decomposition, moment matrices and applications[END_REF]. This approach is closely related to the simultaneous diagonalisation of tensor slices, but follows a more algebraic perspective. Eigenvectors of symmetric tensors have been used to compute their decompositions in [START_REF] Oeding | Eigenvectors of tensors and algorithms for Waring decomposition[END_REF]. In [START_REF] Harmouch | Structured low rank decomposition of multivariate Hankel matrices[END_REF], Singular Value Decomposition and eigenvector computation are used to decompose a symmetric tensor, when its rank is smaller than the smallest size of its Hankel matrix in degree less than half the order of the tensor. In Section 3, we describe a new algorithm, involving Singular Value Decomposition and simultaneous diagonalisation, to compute the decomposition of an identifiable tensor, which interpolation degree is smaller that half the order of the tensor.

Numerical methods such as homotopy continuation have been applied to tensor decomposition in [START_REF] Hauenstein | Homotopy techniques for tensor decomposition and perfect identifiability[END_REF][START_REF] Bernardi | Tensor decomposition and homotopy continuation[END_REF]. Distance minimisation methods to compute low rank approximations of tensors have also been investigated. Alternating Least Squares (ALS) methods, updating alternately the different factors of the tensor decomposition, is a popular approach (see e.g. [CC70, CHLZ12, Har70, KB09]), but suffers from a slow convergence [START_REF] Espig | On the convergence of alternating least squares optimisation in tensor format representations[END_REF][START_REF] Uschmajew | Local convergence of the alternating least squares algorithm for canonical tensor approximation[END_REF]. Other iterative methods such as quasi-Newton methods have been considered to improve the convergence speed. See e.g. [HH82, Paa99, PTC13, SL10, SBL13, TB06, BV18] for multilinear tensors. A Riemannian Newton iteration for symmetric tensors is presented in [START_REF] Khouja | Riemannian newton optimization methods for the symmetric tensor approximation problem[END_REF]. In [START_REF] Kileel | Subspace power method for symmetric tensor decomposition and generalized PCA[END_REF], a method for decomposing real even-order symmetric tensors, called Subspace Power Method (SPM), and similar to the power method for matrix eigenvector computation, is proposed. In these methods, the choice of the initial decomposition is crucial. In the applications of these algorithms, the initial point is often chosen at random, yielding approximate decompositions which can hardly be controlled. Tensor decomposition methods have numerous applications [START_REF] Kolda | Tensor decompositions and applications[END_REF]. Some of them were exploited more recently in Machine Learning. In [START_REF] Hsu | Learning mixtures of spherical gaussians: Moment methods and spectral decompositions[END_REF], symmetric tensor decompositions for moment tensors are studied for spherical Gaussian mixtures. Moment methods have been further investigated for Latent Dirichlet Allocation models, topic or multiview models in [AGH + 14, JGKA19]. In [START_REF] Ruffini | Clustering patients with tensor decomposition[END_REF], a tensor decomposition technique based on Alternate Least Squares (ALS) is used to initialise the Expectation Maximisation (EM) algorithm, for a mixture of discrete distributions (which are not Gaussian distributions).

An overview of tensor decomposition methods in Machine Learning can be found in [START_REF] Rabanser | Introduction to Tensor Decompositions and their Applications in Machine Learning[END_REF].

After reviewing Gaussian mixtures and moment methods in Section 2, we present in Section 3 an algebraic symmetric tensor decomposition method for identifiable tensors. In Section 4, we apply this algorithm for recovering Gaussian mixtures and show its impact on providing good initialisation point in the EM algorithm, in comparison with other state-of-the-art approaches.

Gaussian mixtures and high order moments

In this section, we review Gaussian mixture models and their applications to clustering.

Gaussian mixtures

Suppose that we wish to deal with some Euclidean data x ∈ R m , coming from a population composed of r homogeneous sub-populations (often called clusters). A reasonable assumption is then that each sub-population can be modelled using a simple probability distribution (e.g. Gaussian). This idea is at the heart of the notion of mixture distribution. The prime example of mixture is the Gaussian mixture, whose probability density over R m is defined as

p θ (x) = r j=1 ω j N (x|µ j , Σ j ), (1)
where N (•|µ, Σ) denotes the Gaussian density with mean µ ∈ R m and definite positive covariance matrices Σ ∈ S ++ m . The mixture is parametrised by a typically unknown θ = (ω 1 , ..., ω r , µ 1 , ..., µ r , Σ 1 , ..., Σ r ), composed of

• ω = (ω 1 , ..., ω r ), that belong to the r-simplex and correspond to the cluster proportions,

• µ j and Σ j , that correspond respectively to the mean and covariance of each cluster j ∈ {1, ..., r}.

Gaussian mixtures are ubiquitous objects in statistics and machine learning, and own their popularity to many reasons. Let us briefly mention a few of these.

Density estimation. If r is allowed to be sufficiently large, it is possible to approximate any probability density using a Gaussian mixture (see e.g. [START_REF] Nguyen | Approximation of probability density functions via location-scale finite mixtures in Lebesgue spaces[END_REF]). This motivates the use of Gaussian mixtures as powerful density estimators that can be subsequently used for downstream tasks such as missing data imputation [START_REF] Di Zio | Imputation through finite Gaussian mixture models[END_REF], supervised classification [START_REF] Hastie | Discriminant analysis by Gaussian mixtures[END_REF], or image classification [START_REF] Sánchez | Image classification with the Fisher vector: Theory and practice[END_REF] and denoising [START_REF] Houdard | High-dimensional mixture models for unsupervised image denoising (HDMI)[END_REF].

Clustering. Perhaps the most common use of Gaussian mixtures is clustering, also called unsupervised classification. The task of clustering consists in uncovering homogeneous groups among the data at hand. Within the context of Gaussian mixtures, each group generally corresponds to a single Gaussian distribution, as in Equation (1). If the parameters of a mixture are known, then each point may be clustered using the posterior probabilities obtained via Bayes's rule:

∀x ∈ R m , k ∈ {1, ..., r}, Pr(x belongs to cluster j) = ω j N (x|µ j , Σ j ) p θ (x) . (2)
Detailed reviews on mixture models and their applications, notably to clustering, can be found in [FR02, BCMR19, MLR19].

Learning mixture models

The main statistical question pertaining mixture models is to estimate the parameters θ = (ω 1 , ..., ω r , µ 1 , ..., µ r , Σ 1 , ..., Σ r ) based on a data set x 1 , ..., x n . Typically, X 1 , ..., X n are assumed to be independent and identically distributed random variables with common density p data . The problem of statistical estimation is then to find some θ such that p θ ≈ p data . There are many approaches to this question, the most famous one being the maximum likelihood method. Maximum likelihood is based on the idea that maximising the log-likelihood function

(θ) = n i=1 log p θ (x i ), (3)
will lead to appropriate values of θ. One heuristic reason of the good behaviour of maximum likelihood is that (θ) can be seen as a measure of how likely the observed data is, according to the mixture model p θ . This means that the maximum likelihood estimate will be the value of θ that renders the observed data the likeliest. Another interesting interpretation of maximum likelihood in information-theoretic: when n -→ ∞, maximising the log-likelihood is equivalent to minimising the Kullback-Leibler divergence (an information-theoretic measure of distance between probability distributions) between p θ and p data , thus giving a precise sense to the statement p θ ≈ p data (see e.g. [Bis06, Section 1.6.1]). For more details on the properties of maximum likelihood, see e.g. [VdV98, Section 5.5].

In the specific case of a mixture model, performing maximum-likelihood is however complex for several reasons. Firstly, as shown for instance by [START_REF] Le | Maximum likelihood: an introduction[END_REF], finding a global maximum is actually often ill-posed in the sense that some problematic values of θ will lead to (θ) = ∞ while being very poor models of the data. While focusing on local rather global maxima will fix this first issue in a sense, iterative optimisation algorithms are likely to pursue these unfortunate global maxima. Because of the peculiarities of mixture likelihoods, the most popular algorithm for maximising (θ) is the expectation maximisation (EM, [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]) algorithm, an iterative algorithm specialised for dealing with log-likelihoods of latent variable models. The EM algorithm is usually preferred to more generic gradient-based optimisation algorithms [START_REF] Xu | On convergence properties of the em algorithm for gaussian mixtures[END_REF]. In a nutshell, at each iteration, the EM algorithm clusters the data using Equation (2), and then computes the mean and covariance of each cluster. This iterative scheme is related to another popular clustering algorithm known as k-means (the close relationship between the two algorithms is detailed in [Bis06, Section 9]). A key issue when using the EM algorithm for a Gaussian mixture is the choice of initialisation. Indeed, a poor choice may lead to degenerate solutions, extremely slow convergence, or poor local optima (see [START_REF] Baudry | Em for mixtures[END_REF] and references therein). We will see in this paper that good initial points can be obtained by using another estimation method called the method of moments (as was previously noted by [START_REF] Ruffini | Clustering patients with tensor decomposition[END_REF] in a context of mixtures of multivariate Bernoulli distributions).

The method of moments is a general alternative to maximum likelihood. The idea is to choose several functions g 1 : R m -→ R q 1 , ..., g d : R m -→ R q d called moments, and to find θ by attempting to solve the system of equations

   E x∼p data [g 1 (x)] = E x∼p θ [g 1 (x)] ... E x∼p data [g d (x)] = E x∼p θ [g d (x)]. (4)
Of course, since p data is unknown, solving (4) is not feasible. However, one may replace the expected moments by empirical versions, and solve instead

   1 n n i=1 g 1 (x i ) = E x∼p θ [g 1 (x)] ... 1 n n i=1 g d (x i ) = E x∼p θ [g d (x)].
(5)

A very simple example of this, in the univariate m = 1 case, when g 1 (x) = x, and g 2 (x) = x 2 . Then, solving (4) will ensure that the distributions of the model p θ and the data p data have the same mean and variance. However, many very different distributions have identical mean and variance! A natural refinement of the previous idea is to consider also higher-order moments g 3 (x) = x 3 , g 4 (x) = x 4 , .... This will considerably improve the estimates found using the method of moments. This approach was pioneered by [START_REF] Pearson | Contributions to the mathematical theory of evolution[END_REF] for learning univariate Gaussian mixtures.

In the more general multivariate case m > 1, following [START_REF] Hsu | Learning mixtures of spherical gaussians: Moment methods and spectral decompositions[END_REF], the moments chosen can be tensor products, as we detail in the next section in case of a Gaussian mixture with spherical covariances.

Learning structure from tensor decomposition

In this section, we describe the moment tensors revealing the structure of spherical Gaussian mixtures and how it can be decomposed using standard linear algebra operations.

Let X = (X 1 , . . . , X m ) be a set of variables. The ring of polynomials in

X with coefficients in C is denoted C[X]. The space of homogeneous polynomials of degree d ∈ N is denoted C[X] d .
We recall that a symmetric tensor T of order d (with real coefficients) can be represented by an homogeneous polynomial of degree d in the variables X of the form

T (X) = |α|=d T α d α X α where α = (α 1 , . . . , α n ) ∈ N m , |α| = α 1 + • • • + α m = d, T α ∈ R, d α = d! α 1 !•••αm! , X α = X α 1 1 • • • X αm m .
A decomposition of T as a sum of d th power of linear forms is of the form

T (X) = r i=1 ω i (ξ i • X) d (6)
where

ξ i = (ξ i,1 , . . . , ξ i,m ) ∈ C m and (ξ i • X) = m j=1 ξ i,j X j .
When r is the minimal number of terms in such a decomposition, it is called the rank of T and the decomposition is called a rank decomposition (or a Waring decomposition) of T (X).

We say that the decomposition is unique if the lines spanned by ξ 1 , . . . , ξ r form a unique set of lines with no repetition. In this case, the decomposition of T is unique after normalisation of the vectors ξ i up to permutation (and sign change when d is even). A tensor T with a unique decomposition is called an identifiable tensor. Then the Waring decompositions of T are of the form

T (X) = r i=1 ω i λ -d i (λ i ξ i • X) d for λ i = 0, i ∈ [r]. Given a random variable x ∈ R m , its moments are T α = E[x α 1 1 • • • x αm m ] for α = (α 1 , . . . , α m ) ∈ N m . The symmetric tensor of all moments of order d of x is E[(x • X) d ] = |α|=d E[x α 1 1 • • • x αm m ] d α X α .

The structure of the moment tensor

We aim at recovering the hidden structure a random variable, from the decomposition of its d th order moment tensor. This is possible in some circumstances, that we detail hereafter.

Assumption 3.1. The random variable x ∈ R m is a mixture of spherical Gaussians of probability density (1) with parameters θ = (ω 1 , ..., ω r , µ 1 , ..., µ r , σ 2 1 I m , , ..., σ 2 r I m ) such that r ≤ m. Theorem 3.2 ([HK13]). Under the previous assumption, let

• σ2 be the smallest eigenvalue of E[(x -E[x]) ⊗ (x -E[x])] and v a corresponding unit eigenvector, • M 1 (X) = E[(x • X)(v • (x -E[x])) 2 ], • M 2 (X) = E[(x • X) 2 ] -σ2 X 2 , • M 3 (X) = E[(x • X) 3 ] -3 X 2 M 1 (X). Then σ2 = r i=1 ω i σ 2 i and M 1 (X) = r i=1 ω i σ 2 i (µ i • X), M 2 (X) = r i=1 ω i (µ i • X) 2 , M 3 (X) = r i=1 ω i (µ i • X) 3 . (7)
To analyse the properties of the decomposition (7), we introduce the apolar product on tensors: For two homogeneous polynomials p(X)

= |α|=d d α p α X α and q(X) = |α|=d d α q α X α of degree d, in C[X] d ,
m : ∀u ∈ U m , p(u X), q(u X) d = p(X), q(X) d .
It also satisfies the following properties.

For v ∈ C m , v(X) d = (v • X) d = (v 1 X 1 + • • • + v m X m ) d , p ∈ C[X] d , q ∈ C[X] d-1 , we have : • (v • X) d , p d = p(v), • p, X i q d = 1 d ∂ X i p, q d-1
. For an homogeneous polynomial T of degree d ∈ N (or equivalently a symmetric tensor of order d), we define the Hankel operator of T in degree k ≤ d as the map

H k,d-k T : p ∈ C[X] d-k → [ T, X α p d ] |α|=k ∈ C s k where s k = m+k-1 k = dim C[X] k is the number of monomials of degree k in X. The matrix of H k,d-k T in the basis (X β ) |β|=d-k is H k,d-k T = ( T, X α+β d ) |α|=k,|β|=d-k .
From the properties of the apolar product, we see that H 1,d-1

T : p → 1 d [ ∂ X i T, p d-1 ] 1≤i≤m . For ξ ∈ C m and k ∈ N, let ξ (k) = (ξ α ) |α|=k . We also check that if T = (ξ • X) d with ξ ∈ C m , then H k,d-k (ξ•X) d = ξ (k) ⊗ ξ (d-k
) is of rank 1 and its image is spanned by the vector ξ (k) .

Proposition 3.3. Assume that r ≤ m, w i > 0 for i ∈ [r] and µ 1 , . . . , µ r ∈ R m are linearly independent. The symmetric tensor M 3 (X) is identifiable, of rank r and has a unique Waring decomposition satisfying (7).

Proof. Assume that M 3 (X) has a decomposition of the form (7). Since the vector µ 1 , . . . , µ r are linearly independent, by a linear change of coordinates in Gl m , we can further assume that µ 1 = e 1 , . . . , µ r = e r are the first r vectors of the canonical basis of R m . In this coordinate system, M 3 (X) = r i=1 X 3 i and the matrix H 1,2 M 3 in a convenient basis has a r × r identity block and zero elsewhere. Thus H 1,2 M 3 is of rank r. Its kernel of dimension 1 2 m (m + 1) -r is spanned by the polynomials

X i X j with (i, j) = (k, k) for k ∈ [r]. The kernel of H 1,2
M 3 is thus the space of homogeneous polynomials of degree 2, vanishing at e 1 , . . . , e r ∈ R n .

If M 3 (X) can be decomposed as

M 3 (X) = r i=1 ω i (µ i • X) 3 with ω i ∈ C, µ i ∈ C m and r < r, then H 1,2
M 3 , as a sum of r < r matrices ω i H 1,2 (µ i •X) 3 of rank 1, would be of rank smaller than r < r, which is a contradiction. Thus a minimal decomposition of M 3 (X) is of length r and r is the rank of M 3 (X).

Let us show that the decomposition (7) of M 3 (X) is unique up to a scaling of the vector µ i , i.e. that M 3 (X) is identifiable. For any Waring decomposition M 3 (X) = r i=1 ω i (µ i •X) 3 , the vectors µ 1 , . . . , µ r are linear independant, since

µ i spans imH 1,2 (µ i •X) 3 and H 1,2 M 3 = r i=1 ω i H 1,2 (µ i •X) 3 is of rank r.
As µ 1 , . . . , µ r can be transformed into e 1 , . . . , e r by a linear change of variables, ker H 1,2 M 3 is also the vector space of homogeneous polynomials of degree 2, vanishing at µ 1 , . . . , µ r ∈ C m . Therefore, the set of {µ 1 , . . . , µ r } coincides, up to a scaling, with the set of points {µ 1 , . . . , µ r } of another Waring decomposition of

M 3 (X) = r i=1 ω i (µ i • X) 3 . This shows that M 3 (X) is identifiable.
Therefore, a Waring decomposition of M 3 (X) is of the form M 3 (X) = r i=1 ωi (μ i • X) 3 with ωi = λ -3 ω i , μi = λ i µ i and λ i = 0 for i ∈ [r]. As μ1 , . . . , μr are linearly independent, the homogeneous polynomials (μ 1 • X) 2 , . . . , (μ r • X) 2 are also linearly independant in C[X] 2 (by a linear change of variables, they are equivalent to X 2 1 , . . . , X 2 r ). Consequently, the relation

M 2 (X) = r i=1 ω i (µ i • X) 2 = r i=1 λ i ωi (μ i • X) 2
defines uniquely λ 1 , . . . , λ r , and M 3 (X) has a unique Waring decomposition, which satisfies the relations (7).

Under Assumption 3.1, the hidden structure of the random variable x can thus be recovered using Algorithm 1. This yields the parameters

ω i ∈ R + , µ i ∈ R m , σ i ∈ R + for i ∈ [r] of the Gaussian mixture x.
In the experimentation, the moments involved in the tensors M i will be approximated by empirical moments and we will compute an approximate decomposition of the empirical moment tensor M3 (X).

Decomposition of identifiable tensors

We describe now an important step of the approach, which is computing a Waring decomposition of a tensor. In this section, we consider a tensor T ∈ C[X] d of order d ∈ N with a Waring decomposition of the form T = r i=1 ω i (ξ i • X) d with ω i ∈ C, ξ i ∈ C m , that we recover by linear algebra techniques, under some hypotheses.

Algorithm 1 Recovering the hidden structure of a Gaussian mixture Input: The moment tensors M 1 (X), M 2 (X), M 3 (X).

• Compute a Waring decomposition of M 3 (X) to get ωi ∈ R, μi ∈ R m , i ∈ [r] such that M 3 (X) = r i=1 ωi (μ i • X) 3 . • Solve the system r i=1 ωi (μ i • X) 2 λ i = M 2 (X) to get λ i ∈ R and ω i = λ 3 i ωi ∈ R + , µ i = λ -1 i μi ∈ R m such that M 3 (X) = r i=1 ω i (µ i • X) 3 and M 2 (X) = r i=1 ω i (µ i • X) 2 . • Solve the system r i=1 ω i (µ i • X)σ 2 i = M 1 (X) to get σ 2 i ∈ R + .
Output:

ω i ∈ R + , µ i ∈ R n , σ 2 i ∈ R + for i ∈ [r].
Definition 3.4. The interpolation degree ι(Ξ) of Ξ = {ξ 1 , . . . , ξ r } ⊂ C m is the smallest degree k of a family of homogenous interpolation polynomials u 1 , . . . , u r ∈ C[X] k at the points Ξ (u i (ξ j ) = δ i,j for i, j ∈ [r]).

For any d ≥ ι(Ξ), there exists a family (ũ i ) i∈[r] of interpolation polynomials of degree d, obtained from an interpolation family

(u i ) i∈[r] in degree ι(Ξ) as ũi = (λ•X) d-ι(Ξ) (λ•ξ i ) d-ι(Ξ) u i for a generic λ ∈ C m such that λ • ξ i = 0 for i ∈ [r].
Notice that if the points Ξ = {ξ 1 , . . . , ξ r } are linearly independent (and therefore r ≤ m), then ι(Ξ) = 1 since a family of linear forms interpolating Ξ can be constructed.

If k ≥ ι(Ξ), then the evaluation map e (k)

Ξ : p ∈ C[X] k → (p(ξ 1 ), . . . , p(ξ r )) ∈ C r is surjective.
Its kernel is the space of homogeneous polynomials of degree k vanishing at Ξ. Any supplementary space admits a basis u 1 , . . . , u r , which is an interpolating family for Ξ in degree k. A property of the interpolation degree is the following: Lemma 3.5. For k > ι(Ξ), the common roots of ker e

Ξ = I(Ξ) k is ∪ r i=1 C ξ i .
Hereafter, we show that tensors T such that rank H k,d-k T = r for k > ι(Ξ) + 1 are identifiable and we describe a numerically robust algorithm to compute their Waring decomposition.

Let

U = (U α,j ) |α|=k,j∈[r] ∈ C s k ×r be such that im U = im H k,d-k
T and U i = (U e i +α,j ) |α|=k-1,j∈[r] be the submatrices of U with the rows indexed by the monomials divisible by

X i for i ∈ [m]. Theorem 3.6. Let T ∈ C[X] d with a decomposition T = r i=1 ω i (ξ i • X) d with ω i ∈ C and ξ i = (ξ i,1 , . . . , ξ i,n ) ∈ C m such that rank H k,d-k T = r for some k ∈ [ι(ξ 1 , . . . , ξ r ) + 1, d].
Then T is identifiable of rank r and there exist invertible matrices E ∈ C s k ×s k , F ∈ C r×r such that

E t U i F = ∆ i 0 ( 8
)
with ∆ i = diag( ξ1,i , . . . , ξr,i ) for i ∈ [m]. For any pair (E, F ), which diagonalises simultaneously [U 1 , . . . , U m ] as in (8), there exist unique ω 1 , . . . , ω r ∈ C such that T = r i=1 ω i (ξ i • X) d with ξ i = ((∆ 1 ) i,i , . . . , (∆ m ) i,i ).

Proof. From the decomposition of T , we have for k ≤ d that

H k,d-k T = r i=1 ω i ξ (k) i ⊗ ξ (d-k) i is a linear combination of r Hankel matrices ξ (k) i ⊗ ξ (d-k) i of rank 1.
If T is of rank r < r, then using its decomposition of rank r , H k,d-k T would be of rank ≤ r < r, which is a contradiction. This shows that T is of rank r.

As rank H k,d-k T = r, we deduce that the image of

H k,d-k T is spanned by ξ (k) 1 , . . . , ξ (k)
r and there exists an invertible matrix F ∈ C r×r such that

U F = [ ξ (k) 1 , . . . , ξ(k) r ]
For any polynomial p ∈ C[X] k , which coefficient vector in the monomial basis (X α ) |α|=k is denoted [p], we have [p] t U F = [p( ξ1 ), . . . , p( ξr )] t . This shows that

U ⊥ = {p ∈ C[X] | [p] t U = 0} is ker e (k)
Ξ . By Lemma 3.5 since k ≥ ι( Ξ), the common roots of the homogeneous polynomials in ker e (k) Ξ are the scalar multiples of Ξ. Consequently, the set of lines spanned by the vectors Ξ of a Waring decomposition of T is uniquely determined as the conjugate of the zero locus of

U ⊥ ⊂ C[X] k and T is identifiable. For any p ∈ C[X] k-1 represented by its coefficient vector [p] in the monomial basis (X α ) |α|=k-1 , we have [p] t U i F = [x i p] t U F = [ ξ1,i p( ξ1 ), . . . , ξr,i p( ξr )] t . (9)
Let E be the coefficient matrix of a basis u 1 , . . . , u r , v r+1 , . . . , v s k-1 of C[X] k-1 , such that u 1 , . . . , u r is an interpolating family for Ξ = { ξ1 , . . . , ξr } and v r+1 , . . . , v s k-1 is a basis of ker e (k-1) Ξ .

The matrix E is invertible by construction, and we deduce from (9) that E t U i F = diag( ξ1,i , . . . , ξr,i ) 0 .

Let us show conversely that for any pair of matrices (E , F ), which diagonalises simultaneously [U 1 , . . . , U m ] as in (8) with ∆ i = diag( ξ 1,i , . . . , ξ r,i ), there exist unique ω 1 , . . . , ω r ∈ C such that T

= r i=1 ω i (ξ i • X) d . Let u 1 , . . . , u r , v r+1 , . . . , v s k-1 ∈ C[X] be the polynomials corresponding to the columns of E . Then for a generic λ = (λ 1 , . . . , λ r ) ∈ C m , we have diag((λ • ξ 1 ), . . . , (λ • ξ r )) = m i=1 λ i [u 1 , . . . , u r ] t U i F = m i=1 λ i [u 1 , . . . , u r ] t U i F (F -1 F ) = [(λ • ξj ) u i ( ξj )] i,j∈[r] F -1 F = diag((λ • ξ1 ), . . . , (λ • ξr )) [u i ( ξj )] i,j∈[r] F -1 F .
As λ ∈ C m is generic and λ • ξi = 0 for i ∈ [r], we deduce that ∆ = [u i ( ξj )] i,j∈[r] F -1 F is a diagonal and invertible matrix and that ξ i = ∆i,i ξ i with ∆ i,i = 0. Then we have (

ξ i • X) d = ∆d i,i (ξ i • X) d and T = r i=1 ω i (ξ i • X) d with ω i = ∆-d i,i ω i , which concludes the proof of the theorem.
This leads to Algorithm 2 to compute a Waring decomposition of an identifiable tensor T .

Algorithm 2 Decomposition of an identifiable tensor Input: T ∈ C[X] d , which admits a decomposition with r points Ξ = {ξ 1 , . . . , ξ r } and k > ι(Ξ).

• Compute the Singular Value Decomposition of H

k,d-k T = U S V t ;
• Deduce the rank r of H k,d-k T , take the first r columns of U and build the submatrices U i with rows indexed by the monomials (X i X α ) |α|=k-1 for i ∈ [n];

• Compute a simultaneous diagonalisation of the pencil [U 1 . . . , U m ] as E t U i F = diag( ξ1,i , . . . , ξr,i ) 0 and deduce the points

ξ i = (ξ i,1 , . . . , ξ i,m ) ∈ C m for i ∈ [r];
• Compute the weights ω 1 , . . . , ω r by solving the linear system

T = r i=1 ω i (ξ i • X) d ; Output: ω i ∈ C, ξ i ∈ C m s.t. T = r i=1 ω i (ξ i • X) d .

Numerical experimentations

The model used in this section is the Gaussian Mixture Model (GMM) with differing spherical covariance matrices. Recall that if x = (x 1 , . . . , x n ) is a sample of n independent observations from r multivariate Gaussian mixture with differing spherical covariance matrices of dimension m, and h = (h 1 , h 2 , . . . , h n ) is the latent variable that determine the component from which the observation originates, then:

x i | (h i = k) ∼ N m (µ k , σ 2 k I m ) where, Pr(h i = k) = ω k , for k ∈ [r], such that r k=1 ω k = 1.
The aim of statistical inference is to find the unknown parametrs µ k , σ 2 k and w k , for k ∈ [r] from the data x. This can be done by finding the maximum likelihood estimation (MLE) i.e. finding the optimal maximum of the likelihood function associated to this model. The expectation maximisation algorithm (EM) [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF], usually used for finding MLEs, is an iterative algorithm in which the initialisation i.e. the initial estimation of the latent parameters is crucial, since various initialisations can lead to different local maxima of the likelihood function, consequently, yielding different clustering partition. Thus, in this section we compare the clustering results obtained by different initialisation of the EM algorithm against the initialisation by the method of moments through examples of simulated (subsection 4.1) and real (subsection 4.2) datasets. We fix a maximum of 100 iterations of the EM algorithm. The different initialisation considered in this section are the following:

• The k-means method [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] according to the following strategy:

The best partition obtained out of 50 runs of the k-means algorithm.

• The method of moments, where Algorithm 1 is applied to build the moments and Algorithm 2 is applied to the empirical moment tensor corresponding to M 3 (X) (see Theorem 3.2), with less than 5 Riemannian Newton iterations [START_REF] Khouja | Riemannian newton optimization methods for the symmetric tensor approximation problem[END_REF] to reduce the distance between the empirical moment tensor and its decomposition.

• The Model-based hierarchical agglomerative clustering algorithm (MBHC) [START_REF] Vaithyanathan | Model-based hierarchical clustering[END_REF][START_REF] Fraley | Algorithms for model-based gaussian hierarchical clustering[END_REF].

• The emEM strategy [START_REF] Biernacki | Choosing starting values for the em algorithm for getting the highest likelihood in multivariate gaussian mixture models[END_REF] as in [LIL + 15] which makes 5 iterations for each of 50 short runs of EM, and follows the one which maximises the log-likelihood function by a long run of EM.

The k-means, MBHC and emEM are common strategies for initialising the EM algorithm for GMMs. The comparison among the different EM initialisation strategies is based on three measures: The Bayesian Information Criterion (BIC) [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF][START_REF] Fraley | How many clusters? which clustering method? answers via model-based cluster analysis[END_REF], the Adjusted Rand Index (ARI) [START_REF] Hubert | Comparing partitions[END_REF], and the error rate (errorRate). The BIC is a penalized-likelihood criterion given by the following formula

BIC = -2 ( θ) + log(n)ν,
where is the log-likelihood function , θ is the MLE which maximises the log-likelihood function and ν is the number of the estimated parameters. This criterion measures the quality of the model such that for comparing models the one with the largest BIC value among the other models is the most fitted to the studied dataset. The ARI criterion measures the similarity between the estimated clustering obtained by the applied model and the exact true clustering. Its value is bounded between 0 and 1. The more this measure is close to 1 the more the estimated clustering is accurate. The error rate measure can be viewed as an alternative of the ARI. In fact this criterion measures the minimum error between the predicted clustering and the true clustering, and thus low error rate means high agreement between the estimated and the true clustering. The former criteria as well as the EM algorithm are used from the tools of the package mclust [START_REF] Scrucca | mclust 5: clustering, classification and density estimation using Gaussian finite mixture models[END_REF] in R programming language.

Simulation

We performed 100 simulations from each of the two models described in examples 4.1 and 4.2. We counted the instances where each of the considered initialising strategies for the EM could find throughout the 100 simulated data and among the other initialisation methods the largest BIC, the highest ARI, ARI≥ 0.99 (as in this case the clustering obtained is the most accurate) and the lowest errorRate. The values of the BIC, ARI, errorRate and consumed time of the different considered initialisation strategies for one dataset sampled according to the model of Example 4.1 (resp. 4.2) are presented in Table 1 (resp. 3), and Figure 1 (resp. 2) shows a two-dimensional visualisation of the observations according to the first four features, the observations in the upper panels are labeled according to the actual clustering, while they are labeled in the lower panels according to the clustering obtained by the EM algorithm initialised by the method of moments. In order to have an estimation about the numerical stability of the obtained results, we repeat the same numerical experiment for each example 20 times and we compute the means (Table 2, 4) and the variances (values in parentheses in Table 2, 4) of the 20 percentages obtained of each of the BIC, ARI, ARI≥ 0.99 and errorRate values for the different initialising strategies. As we mentioned before the initialisation strategies considered in this comparison against the method of moments are common and have, in general, good numerical behavior. Nevertheless, we cannot expect all the initialisation strategies that exist for the EM algorithm to work well in all the cases [START_REF] Biernacki | Choosing starting values for the em algorithm for getting the highest likelihood in multivariate gaussian mixture models[END_REF][START_REF] Melnykov | Finite mixture models and model-based clustering[END_REF]. Hereafter, two examples are chosen in such a way to present some cases where the common initialisation strategies k-means, MBHC and emEM have some difficulties to provide a good initialisation to the EM algorithm for the GMMs with differing spherical covariance matrices, or in other words where the initialisation by the method of moments outperforms the other considered initialisations. For instance, we put in each of these two examples one cluster of small size (the blue cluster in Figure 1, the red cluster in Figure 2), we want to make the clusters overlap, since these initialisation strategies could misscluster the dataset if the clusters are intersecting. We notice that this choice of the mean vectors and the different variances in each of the two examples yields a dataset with the expected clustering characteristic.

Example 4.1. In the first simulation example, a multivariate dataset (m=6) of n=1000 observations generated with r=4 clusters according to the following parameters:

• The probability vector: ω = (0.2782, 0.0139, 0.3324, 0.3756) T .

• The mean vectors: µ 1 = (-5.0, -9.0, 8.0, 8.0, 2.0, 5.0) T , µ 2 = (-7.0, 6.0, -1.0, 6.0, -8.0, -10.0) T , µ 3 = (-4.0, -10.0, -5.0, 1.0, 5.0, 4.0) T , µ 4 = (-6.0, 6.0, 5.0, 4.0, -1.0, -1.0) T .

• The variances: σ 2 1 = 1.5, σ 2 2 = 2.5, σ 2 3 = 5.0, σ 2 4 = 15.0. Example 4.2. In the second simulation example, a multivariate dataset (m=5) of n=1000 observations generated with r=3 clusters according to the following parameters:

• The probability vector: ω = (0.0930, 0.2151, 0.6918) T .

• The mean vectors: µ 1 = (7.0, -4.0, -4.0, -6.0, -4.0) T , µ 2 = (2.0, -4.0, -6.0, -10.0, -3.0) T , µ 3 = (4.0, -4.0, -5.0, 6.0, 1.0) T .

• The variances: σ 2 1 = 5.0, σ 2 2 = 10.0, σ 2 3 = 15.0. The Table 2, 4 show that in Example 4.1, 4.2 the best results among the considered initialising strategies are for the method of moments. In fact, in the former two tables we see that the method of moments found throughout the 100 simulated datasets, in average (by runing the numerical experiment 20 times), the largest BIC, highest ARI, ARI≥ 0.99 and lowest errorRate among the other initialisation strategies in more instances than all the other considered initialisation method, implying in this context marked outperformance for the moments initialisation method. Note that the consumed time (see. Table 1, 3) tends to be higher in the method of moments than in the other initialisation strategies. This is expected since stochastic approaches (to which the methods k-means, MBHC and emEM belong) outperform the deterministic approaches (as the method of moments) in this term.

Real data

In this subsection we present four examples of real datasets, for which we know already their number of clusters, and we report the different BIC, ARI and errorRate values as well as the consumed time attained by the EM algorithm initialised by the different considered initialisation strategies and used with the GMM of different spherical covariance matrices. The explored real data are: The famous iris data [Fis36, DT17] widely used as an example of clustering to test the algorithms, Diabetes [START_REF] Reaven | An attempt to define the nature of chemical diabetes using a multidimensional analysis[END_REF], olive oil [START_REF] Azzalini | Clustering via nonparametric density estimation: The R package pdfcluster[END_REF], and MNIST [START_REF] Deng | The mnist database of handwritten digit images for machine learning research[END_REF].

Example 4.3 (Iris). The iris dataset contains four physical measurements (length and width of sepals and petals) for 50 samples of three species of iris (setosa, virginica and versicolor). The number of features is m = 4 and the number of clusters is r = 3. The four initialisation strategies yield the same BIC value. The ARI and the errorRate values Example 4.4 (Diabetes). The Diabete dataset [RM79] contains three measurements: glucose, insulin and sspg; made on 145 non-obese adult patients classified into three types of diabetes: Normal, Overt, and Chemical. Herein, in this example m = r = 3. We apply the different initialisation strategies for the EM algorithm, the Table 6 shows the results. Despite the fact that k-means method is the fastest method in this example, the ARI and the BIC are noticeably lower than in the other methods. Concerning the method of moments, it succeeds to have quite similar scores to the other methods in this example, but with a bigger computation time.

Example 4.5 (Olive oil). The olive oil data set contains the chemical composition (8 chemical properties) of 572 olive oils. They are derived from three different macro-areas in Italy (South, Sardinia and Centre North). The dataset contains nine regions from which the olive oils were taken in Italy. Thus we can cluster this dataset according to the macro-areas (r = 3) or the region (r = 9). As the number of features in this dataset is m = 8, we choose r = 3, so that the This shows that for these datasets which are not well fitted by the mixture of spherical Gaussians, the moment method can still give good initialisations for the EM algorithm, in comparision with the common initialisation strategies.

Example 4.6 (MNIST digit image database). The MNIST digit image database [START_REF] Deng | The mnist database of handwritten digit images for machine learning research[END_REF] is a large database that contains images of 28 × 28 pixels for handwritten digits (0 to 9). Each pixel contains an integer between 0 and 255 that represents the grayscale levels. The number of features is 28 × 28 = 784. We choose the MNIST digit image dataset which contains 60000 images. We take a subset of this dataset that contains the images of label 0 or 1. The size of the subset is 12665 images. Since the number of features is quite large (784), and we aim to test a spherical Gaussian mixture model, a good practice in this case is to apply one of the dimensionality reduction strategies. Roughly speaking, the dimensionality reduction strategies aim to reduce the number of features such that a high percentage of the information within the dataset is conserved. In other words, the performance in term of accuracy of the clustering methods will not be noticeably affected by this reduction, and on the other hand this will reduce considerably the time of computation. For this purpose, we choose to apply the Principal Component Analysis transformation (PCA) [F.R01, Jol11]. We conserve the first five variables given by this transformation (see Figure 3). The dataset that we consider in this example contains 12665 observations, the number of clusters is r = 2, and the number of features is m = 5. We apply the different initialisation strategies and we report the results in Table 8. As we can see, the results given by the method of moments in Table 8 are very satisfactory in comparison with the other initialisation strategies with ARI= 0.9308. In particular, the method of moments clearly outperforms MBHC method in this regard, in term of accuracy and the time of computation. In fact, the MBHC takes 543.4 seconds without reaching a good ARI score. This example sheds some light on the performance of the method of moments. The large number of samples (in this example equal to 12665) does not have a high impact on the computation time, which is not the case, for the MBHC method, where this factor increases significantly its computation time. Moreover, it is true that a large number of features could have a negative impact on the computation time of the method of moments, but it is not a sever limitation since as we saw in this example, this can be efficiently remedied by applying one of the dimensionality reduction techniques. In this regard, some recent work [START_REF] João | Tensor moments of gaussian mixture models: Theory and applications[END_REF] studies how the computation complexity of the moment method can be reduced while conserving its desirable high accuracy property. Conducting more research in this direction, we believe that the method of moments will have more sophisticated and competitive (in term of computation time) developments in the future.

Conclusion

In the context of unsupervised machine learning, the type of models to be recovered plays an important role. For Gaussian mixture models, where iterative methods such as Expectation Maximisation algorithms are applied, the choice of the initialisation is also crucial to recover an accurate model of a given dataset. We demonstrated in the experimentation that tensor decomposition techniques can provide a good initial point for the EM algorithm, and that the moment tensor method outperforms the other state-of-the-art strategies, when datasets are well represented by spherical Gaussian mixture models. For that purpose, we presented a new tensor decomposition algorithm adapted to the decomposition of identifiable tensors with low interpolation degree, which applies to a 3 rd order moment tensors associated to the data distribution as we have shown.

  their apolar product is p, q d := |α|=d d α pα q α . The apolar norm of p is ||p|| d = p, p d = |α|=d d α pα p α . The apolar product is invariant by a linear change of variables of the unitary group U

Ξ

  is the union ∪ r i=1 C ξ i of lines spanned by ξ 1 , . . . , ξ r ∈ C m . Proof. As ι(Ξ) + 1 is the Castelnuovo-Mumford regularity of the vanishing ideal I(Ξ) = {p ∈ C[X] | p homogeneous, p(ξ) = 0 for ξ ∈ Ξ} [Eis05][Ch.4], it is generated in degree k > ι(Ξ) and the common roots of ker e (k)

Figure 1 :

 1 Figure 1: Scatterplot matrix for the sampled dataset of Example 4.1 projected onto the first four variables (features): upper panels show scatterplots for pairs of variables in the original clustering; lower panels show the clustering obtained by applying the EM algorithm initialised by the method of moments.

Figure 2 :

 2 Figure 2: Scatterplot matrix for the sampled dataset of Example 4.2 projected onto the first four variables (features): upper panels show scatterplots for pairs of variables in the original clustering; lower panels show the clustering obtained by applying the EM algorithm initialised by the method of moments.

Figure 3 :

 3 Figure 3: Scatterplot for pairs of variables: upper panels show the first five features obtained by applying the PCA transformation on the dataset of Example 4.6. The graphs points marked according to the true two classes 0 and 1.

Table 1 :

 1 Numerical results of one data set of Example 4.1

	Method	BIC	ARI errorRate time(s)
	em km	-29590.48 0.8281	0.168	0.045
	em mom -29492.11	1.0	0.0	0.547
	em mbhc	-29594.97 0.8574	0.099	0.287
	em emEM -29593.18 0.8366	0.132	0.171

Table 2 :

 2 Estimation of the stability of Example 4.1 results

	Method	BIC	ARI	ARI ≥ 0.99	errorRate
	em km	38.35% (37.82) 47.6% (21.41) 48.85% (21.61)	47.6% (21.2)
	em mom 74.8% (41.01) 88.75% (15.36) 83.4% (18.36) 88.60% (14.46)
	em mbhc 10.75% (12.41) 15.9% (17.57) 15.55% (22.99) 15.9% (19.46)
	em emEM	7.3% (8.43)	14.5% (8.05)	12.6% (17.83)	14.95% (7.52)

Table 3 :

 3 Numerical results of one data set of Example 4.2

	Method	BIC	ARI	errorRate time(s)
	em km	-28360.30	0.4352	0.309	0.051
	em mom -28246.02 0.9498	0.03	0.504
	em mbhc	-28358.67	0.3197	0.384	0.292
	em emEM -28360.42	0.4408	0.296	0.141

Table 4 :

 4 Estimation of the stability of Example 4.2 results

	Method	BIC	ARI	ARI ≥ 0.99	errorRate
	em km	0.45% (0.576)	0.05% (0.05)	0.0% (0.0)	0.1%(0.095)
	em mom 50.0% (18.63) 92.35% (9.82) 0.0% (0.0) 92.1% (7.46)
	em mbhc 49.35% (19.82) 2.45% (3.63)	0.0% (0.0) 2.45% (2.58)
	em emEM 0.3% (0.326)	5.2% (4.48)	0.0% (0.0)	5.9% (5.36)

Table 5 :

 5 Numerical results of Example 4.3

	Method	BIC	ARI	errorRate time(s)
	em km	-1227.6656	0.6199	0.167	0.007
	em mom -1227.6676 0.6410	0.153	0.203
	em mbhc	-1227.6696	0.6199	0.167	0.007
	em emEM -1227.6495	0.6302	0.160	0.045

are slightly better with the moment initialisation among the other considered initialisation strategies. On the other hand, the consumed time is clear higher in the moment method initialisation.

Table 6 :

 6 Numerical results of Example 4.4

	Method	BIC	ARI	errorRate time(s)
	em km	-5363.06 0.3371	0.289	0.007
	em mom	-5222.11 0.6355	0.144	0.380
	em mbhc -5221.32 0.6355	0.144	0.008
	em emEM -5221.33 0.6207	0.151	0.049

Table 7 :

 7 Numerical results of Example 4.5 ≤ m for the method of moment is verified. The results show that the MBHC initialisation strategy yields the largest BIC, the highest ARI and the lowest errorRate values among the other initialisation strategies. Nevertheless, the initialisation by the moment method comes in second position after the MBHC strategy in terms of the BIC, ARI and errorRate values, while the K-means and the emEM initialisation strategies attain almost the same values of the previously mentionned criteria.

	Method	BIC	ARI	errorRate time(s)
	em km	-10948.64	0.4018	0.262	0.021
	em mom	-10946.46	0.4532	0.210	0.508
	em mbhc -10625.59 0.5003	0.185	0.080
	em emEM -10948.72	0.4040	0.260	0.087
	condition r

Table 8 :

 8 Numerical results of Example 4.5

	Method	BIC	ARI	errorRate time(s)
	em km	-384977.3	0.9304	0.017	0.537
	em mom	-384978.2 0.9308	0.017	1.87
	em mbhc -382746.2 0.2445	0.252	543.4
	em emEM -384977.6	0.9301 0.0177655	1.80

Acknowledgement. We would like to thank the anonymous reviewers for their valuable remarks that helped us improving the paper.