Revisiting clustering as matrix factorisation on the Stiefel manifold

TL;DR

This paper introduces a novel Bayesian clustering method based on matrix factorization on the Stiefel manifold, providing theoretical guarantees and an efficient sampling algorithm for high-dimensional data.

Contribution

It reformulates clustering as low rank matrix estimation on the Stiefel manifold using Burer-Monteiro factorization, with new prediction bounds and a Langevin sampler.

Findings

Proposes a new Bayesian estimator with theoretical prediction bounds.

Develops a Langevin sampler for efficient computation.

Demonstrates effectiveness on high-dimensional data.

Abstract

This paper studies clustering for possibly high dimensional data (e.g. images, time series, gene expression data, and many other settings), and rephrase it as low rank matrix estimation in the PAC-Bayesian framework. Our approach leverages the well known Burer-Monteiro factorisation strategy from large scale optimisation, in the context of low rank estimation. Moreover, our Burer-Monteiro factors are shown to lie on a Stiefel manifold. We propose a new generalized Bayesian estimator for this problem and prove novel prediction bounds for clustering. We also devise a componentwise Langevin sampler on the Stiefel manifold to compute this estimator.

Equations

$$X=\left[x_1,\dots ,x_n\right].$$

$$M=\left[{\mu }_{k_1},\dots ,{\mu }_{k_n}\right],$$

$$E=\left[{ϵ}_1,\dots ,{ϵ}_n\right]$$

$$T_{i,j}^{\ast }$$

$$T^{\ast }=U^{\ast }{U^{\ast }}^t$$

$$U^{\ast }=\left[\frac{1}{\sqrt{n_1}}{1}_{I_1},\dots ,\frac{1}{\sqrt{n_K}}{1}_{I_K}\right],$$

$$U$$

$${\pi }_{U,O}(U,O)$$

$${\pi }_{U_{i,r}\mid O}(U)$$

$$L_{\lambda }(U)$$

$$\rho (U,O)$$

$${\hat{U}}_{\lambda }$$

$${\nu }_{\min }\le \underset{\tilde{U}\in {\mathbb{O}}_{n,R,+},M\tilde{U}{\tilde{U}}^t=M}{\min }\Vert E(I-\tilde{U}{\tilde{U}}^t){\Vert }_F$$

$${\nu }_{\max }\ge \underset{\tilde{U}\in {\mathbb{O}}_{n,R,+},M\tilde{U}{\tilde{U}}^t=M}{\max }\Vert E(I-\tilde{U}{\tilde{U}}^t){\Vert }_F.$$

$$c_U\left(2+c_U\right)\le ϵ\frac{{\nu }_{\min }}{\Vert M\Vert +\Vert E\Vert },$$

$$\mathbb{E}\left[\left\|M\left(T^{*}-\hat{U}_{\lambda}\hat{U}_{\lambda}^{t}\right)\right\|_{F}\right]\leq(1+\epsilon)\min_{\tilde{U}\in\mathbb{O}_{n,R,+},\ M\tilde{U}\tilde{U}^{t}=M}\quad\Bigg{\{}\|M(T^{*}-\tilde{U}\tilde{U}^{t})\|_{F}\\ +\sqrt{\frac{1}{\exp\left(\frac{\left(\sqrt{c_{U}^{2}-c_{O}^{2}}-\sqrt{nR}\right)^{2}}{2}\right)-1}}\\ +\sqrt{\left(nR-\frac{1}{2}(R^{2}+R)\right)\log(\rho^{-1})+\log(c^{-1})}+(2+\epsilon)\nu_{\max}\Bigg{\}}.$$

$$c_U\left(2+c_U\right)\le ϵ\frac{{\nu }_{\min }}{\sqrt{\left({\max }_{k=1}^K{n}_k\right)\mu (K-1)+1}+{\sigma }_{\max }\left(\sqrt{n}+2\sqrt{d}\right)},$$

$$1-\exp (-d)-\exp (-n{u}^2/8)-K,$$

$$K$$

$$\begin{gathered} \sum _{k=1}^K\sum _{i_k\in I_k}{\left(\sum _{i_k^{\prime }\in I_k}{T}_{\pi ,i_k^{\prime },i_k}^{\ast }-{\hat{U}}_{\lambda ,\pi ,i_k^{\prime }}{\hat{U}}_{\lambda ,\pi ,i_k}^t\right)}^2 \\ \le \frac{(1+ϵ)\sqrt{1+\mu (K-1)}}{1-\mu (K-1)}\underset{\tilde{U}\in {\mathbb{O}}_{n,R,+},M\tilde{U}{\tilde{U}}^t=M}{\min }\Vert M(T^{\ast }-\tilde{U}{\tilde{U}}^t){\Vert }_F \\ +\sqrt{\frac{1}{\exp \left(\frac{{\left(\sqrt{c_U^2-c_O^2}-\sqrt{dR}\right)}^2}{2}\right)-1}}+\sqrt{(dR-\frac{1}{2}(R^2+R))\log ({\rho }^{-1})+\log (c^{-1})} \\ +(2+ϵ){\nu }_{\max }. \end{gathered}$$

$$t_{\min }$$

$$t_{\max }={\left(\frac{{\nu }_{\max }}{{\sigma }_{\max }}-4ϵ\sqrt{nd+u}\right)}^2-\sqrt{n(n-R)},$$

$$[O]$$

$$T_{O}(\mathbb{O}_{d,R})=\Bigg{\{}O\begin{bmatrix}A&-B^{t}\\ B&0\end{bmatrix},\quad\mathrm{with}\quad A\in\mathbb{R}^{R\times R}\ \mathrm{skew\ symmetric}\Bigg{\}}.$$

$$g_c$$

$$\nabla f(O)$$

$$\nabla {\left(\Vert X-X\cdot {\cdot }^t{\Vert }_F^2\right)}_U=\left((X-XU{U}^t{)}^{t}X+X^t(X-XU{U}^t)\right)U.$$

$$O^{(\ell +1)}$$

$$-O^{(\ell )}{\left(\mathrm{sign}(O^{(\ell )})\odot \left(U^{(\ell )}-|O^{(\ell )}|\right)\right)}^t{O}^{(\ell )}$$

$$\displaystyle+\sqrt{2\ h}\ Z_{O}^{(\ell)}\Big{)}\Bigg{)}$$

$$U^{(\ell +1)}$$

$$\displaystyle+X^{t}(X-XU^{(\ell)}U^{(\ell)^{t}})\bigg{)}U^{(\ell)}$$

Full text

Keywords. Clustering, concentration inequalities, non-negative matrix factorisation, Gaussian mixtures, PAC-Bayes, optimisation on manifolds.

1 Introduction

Clustering, i.e., unsupervised classification, is a central problem in machine learning and has attracted great attention since the origins of statistics, via model-based learning, but recently regained a lot of interest from theoreticians, due to its similarities with community detection Arias-Castro and Verzelen [2014], Verzelen and Arias-Castro [2015]. On the application side, clustering is pervasive in data science, and has become a basic tool in computer science, bio-informatics, finance, metrology, to name but a few.

1.1 Historical background

The problem of identifying clusters in a data set can be addressed using an wide variety of tools. Two main approaches can be delineated, namely the model-based approach and the non-model based approach. Techniques such as hierarchical clustering Hastie et al. [2009], minimum spanning tree-based approaches Blum et al. [2016], $K$-means algorithms Hastie et al. [2009], belong to the non-model based family of methods. Model-based techniques mostly rely on mixture modelling McLachlan and Peel [2004] and often offer better interpretability whilst being easily amenable to uncertainty quantification analysis. The EM algorithm Dempster et al. [1977], McLachlan and Peel [2004] is often the algorithm of choice in the frequentist approach while many Monte Carlo Markov Chain techniques have been proposed for estimation in the Bayesian setting.

In recent years, the clustering problem has revived a surge of interest in a different setting, namely community detection in random graphs. Tools from spectral graph theory and convex optimisation, combined with recent breakthrough from random matrix theory were put to work in devising efficient clustering methods that operate in polynomial time. The celebrated example of Max-Cut, a well known NP-hard combinatorial optimisation problem strongly related to bi-clustering and with a tight Semi-Definite Programming (SDP) relaxation discovered by Goemans and Williamson [1995], is an example among the many successes of the convex optimisation approach to addressing machine learning problems. SDP is the class of optimisation problems that consist in minimising a linear function over the sets of Positive Semi-Definite matrices that satisfy a set of linear (in)equalities. Goemans and Williamson [1995] subsequently triggered a long lasting trend of research in convex relaxation with many application in data science, and recent results proposing tighter relaxations to the clustering problem can be found in Guédon and Vershynin [2016], Chrétien et al. [2016], Giraud and Verzelen [2018]. Some of these methods even apply to any kind of data set endowed with a relevant affinity measure computed from the pairwise distances between the data points, and share the common feature of using low-rank matrix representations of the clustering problem. The theoretical tools behind the analysing of the performance of these convex optimisation-based methods are also quite fascinating and range from random matrix theory Bandeira [2018], Vershynin [2018], concentration inequalities for quadratic forms of random vectors Rudelson and Vershynin [2013] and optimisation theory [optimality conditions, see Royer, 2017], localisation arguments in statistical learning theory Giraud and Verzelen [2018], Grothendieck’s inequality Guédon and Vershynin [2016], Montanari and Sen [2015], to name but a few.

The main drawback, however, of the current lines of approach to the performance analysis of these powerful convex SDP and spectral relaxations is that they all depend on the separation between clusters, i.e., the minimum distance between two points from different clusters, a crucial parameter in the aforecited analyses. In real data sets however, sufficient inter-cluster separation rarely holds and overlaps between clusters are the common situation. This leaves open the difficult problem of finding an alternative theoretical route for controlling the estimation error. On the computational side, the sample size is also a problem for SDP relaxations for which off-the-shelf software does not scale to big data. A remedy to this problem is to use the Burer-Monteiro factorisation consisting in solving in $U$ where $X=UU^{t}$ is the variable of the SDP at hand Burer and Monteiro [2003]. The Burer-Monteiro factorisation results in a non-convex optimisation problem whose local minimisers are global minimisers when the number of columns of $U$ is sufficiently large Burer and Monteiro [2005], Boumal et al. [2016]. In practice however, the rank of the sought matrix is simply equal to the number of clusters, and whether such small priors on the rank of the Burer-Monteiro factorisation are compatible with the local/global equivalence of the minimisers in general remains an open question to this day. A final source of frustration in our list, is that there does not seem to exist any method for quantifying the uncertainty of the results in these convex optimisation-based approaches to clustering.

In the present paper, we propose a generalized Bayesian approach to clustering which hinges on low rank estimation of a clustering matrix. We then leverage arguments from the PAC-Bayesian theory for controlling the error which does not use any prior estimate of separation. Our approach is based on the estimation of a normalised version $T^{*}$ of the adjacency matrix of the clustering, which can be factorised into $T^{*}=U^{*}U^{*^{t}}$, where $U^{*}$ has orthonormal, non-negative columns. Leveraging this structure leads to sampling on the intersection of the Stiefel manifold Edelman et al. [1998] and the non-negative orthant, which is another surprising manifestation of the power of non-negative matrix factorisation (NMF) in clustering problems. Solving this factorised version in the PAC-Bayesian setting is the sampling counterpart of the Burer-Monteiro approach to the numerical solution of high dimensional SDP. The PAC-Bayesian approach (initiated by Shawe-Taylor and Williamson [1997], McAllester [1998, 1999], Catoni [2004, 2007]; see Guedj [2019] for a recent survey) moreover makes no prior use of the separation and at the same time makes it possible to obtain state-of-the-art risk bounds.

1.2 Our contribution

The main goal of the present paper is to study the clustering problem from a low rank Stiefel matrix, i.e. matrices with orthonormal columns, view point, and present a PAC-Bayesian analysis of the related statistical estimation problem. Our approach is in particular inspired by recent work on low rank approximation for $k$-means Boutsidis et al. [2009], Cohen et al. [2015], where the representation of clustering using the matrix $T^{*}$ is explicitly stated (although no algorithm is provided), and PAC-Bayesian bounds for Non-Negative Matrix factorisation [as introduced by Alquier and Guedj, 2017, although they do not establish the link between NMF and clustering]. To the best of our knowledge, the representation in Boutsidis et al. [2009] using the matrix $T^{*}$ has never been studied from a statistical learning perspective.

We present our main result (Theorem 1, which states an inequality holding in expectation on the prediction performance) in Section 2. Our second main result is Theorem 2, which specifies the results of Theorem 1 in the case where we assume that the family of means is incoherent. Section 3 is devoted to our algorithm (an alternating Langevin sampler which relies on computing gradients on the Stiefel manifold).

Our choice of the Langevin sampler is motivated by its faster convergence than Metropolis type samplers due to the absence of any rejection step. This is even more true in the case where one needs to account for the manifold structure of the model as is the case in the present paper.

Additional proofs are gathered in the supplementary material.

1.3 Notation

The notation used in the present paper is fairly standard. The canonical scalar product in $\mathbb{R}^{d}$ will be denoted by $\langle\cdot,\cdot\rangle$, the $\ell_{p}$ norms by $\|\cdots\|_{p}$. For matrices in $\mathbb{R}^{d\times n}$, the operator norm will be denoted by $\|\cdot\|$ and the Frobenius norm by $\|\cdot\|_{F}$. The Stiefel manifold of order $(n,R)$, i.e. the set of matrices in $\mathbb{R}^{n\times R}$ with orthonormal columns, will be denoted by $\mathbb{O}_{n,R}$, and $\mathbb{O}_{n,R,+}$ will denote the subset of the Stiefel manifold $\mathbb{O}_{n,R}$ consisting of componentwise nonnegative matrices. The matrices in $\mathbb{O}_{n,R}$ will sometimes be identified with matrices in $\mathbb{R}^{n\times n}$ where the first $R$ columns form an orthonormal family and the remaining $n-R$ columns are set to zero. The gradient operator acting on differentiable multivariate functions will be denoted by $\nabla$.

2 Non-negative factorisation of the Stiefel manifold

This section is devoted to the presentation of our framework and our main theoretical result.

2.1 Model

Let data points $x_{1},\ldots,x_{n}$ be vectors in $\mathbb{R}^{d}$ and let $X$ denote the matrix

[TABLE]

Let $\mu_{1},\ldots,\mu_{K}$ be $K\in\mathbb{N}\backslash\{0\}$ vectors in $\mathbb{R}^{d}$. We will say that $x_{i}$ belongs to cluster $k\in\{1,\dots,K\}$ if $x_{i}=\mu_{k}+E_{i}$ for some centered random vector $E_{i}\in\mathbb{R}^{d}$. For each $i=1,\ldots,n$, we will denote by $k_{i}$ the label of the cluster to which $x_{i}$ belongs. For each $k$, we will denote by $I_{k}$ the index set of the points which belong to cluster $k$ and $n_{k}$ its cardinality. Now, we can decompose $X$ as $X=M+E$ with

[TABLE]

[TABLE]

More assumptions about the noise matrix $E$ will be introduced and their consequences on the performance of our clustering method will be studied in Theorem 2.

2.2 Ideal solution

If we try to estimate the columns of $M$, one simple way is to use a convex combination of the $x_{i}$’s for each of them. In other words, one might try to approximate $X$ by $XT^{*}$ where $T^{*}$ is a $\mathbb{R}^{n\times n}$ matrix. One simple way to proceed is to set $T$ as the matrix which computes the cluster means, given by

[TABLE]

Thus, each column $i=1,\ldots,n$ of $XT^{*}$ is simply the mean over cluster $k_{i}$. This type of solution is well-motivated by the fact that the mean is the least-squares solution of the approximation of $\mu_{k}$ by the observation points. The matrix $T^{*}$ defined as above enjoys the following desirable properties: (i) its rank is exactly the number of clusters (ii) it is nonnegative (iii) the columns corresponding to different clusters are orthogonal.

One important fact to notice is that the eigenvalue decomposition of $T^{*}$ is explicit and given by

[TABLE]

with

[TABLE]

and therefore, all the eigenvalues of $T^{*}$ are equal to one.

Based on this decomposition, we can now focus on estimating $U^{*}$ rather than $T^{*}$, the reason being that working on estimating $U^{*}$ with $\hat{U}\geq 0$ will automatically enforce positive semi-definiteness of $\hat{T}$ (the estimator of $T^{*}$) and non-negativity of its components. Moreover, enforcing the orthogonality of the columns of $\hat{U}$, combined with the non-negativity of its components, will enforce the columns of $\hat{U}$ to have disjoint supports.

Adopting a generalized Bayesian strategy [inspired by Alquier and Guedj, 2017], we will then define a prior distribution on $\hat{U}$ and study the main properties of the resulting (generalized) posterior distribution.

2.3 The latent variable model

In order to perform an accurate estimation, we need to devise meaningful priors which will account for the main constraints our estimator should satisfy, namely (i) nonnegativity of the entries (ii) orthogonality of the columns (iii) the columns have unit $\ell_{2}$ norm (iv) group sparsity of the columns.

In order to simplify this task, we will introduce a latent (matrix) variable $O$ with uniform distribution on the orthogonal group, and build priors on $U$ that will promote group sparsity of the columns and non-negativity (component-wise).

Prior on $(U,O)$.

Let $\mathcal{U}_{R}$ denote the set of matrices of the form

[TABLE]

Let $\mathbb{O}_{n,R}$ denote the Stiefel manifold, i.e., the manifold of all matrices with $R$ orthonormal columns in $\mathbb{R}^{n}$. The prior on $(U,O)\in\mathcal{U}_{R}\times\mathbb{O}_{n,R}$ is given by

[TABLE]

with

[TABLE]

for $i=1,\ldots,n$, and $r=1,\ldots,R$, with $R$ being a fixed integer and $\pi_{O}$ being the uniform distribution on the Stiefel manifold $\mathbb{O}_{n,R}$.

2.4 Generalized posterior and estimator

Following the approach of Alquier and Guedj [2017], we use a loss term [instead of a likelihood, hence the term ”generalized Bayes”, see Guedj, 2019, for a survey] given by

[TABLE]

for some fixed positive parameters $\lambda$ and $\mu$. The resulting generalized posterior (also known as a Gibbs measure) is defined as

[TABLE]

where $Z_{\lambda}$ denotes the normalisation constant $Z_{\lambda}=\int\ L_{\lambda}(U)\ \pi_{U\mid O}(U)\ \pi_{O}(O)\ dU.$ Finally we let $\hat{U}_{\lambda}$ denote the posterior mean of $U$, i.e.

[TABLE]

2.5 A PAC-Bayesian-flavored error bound

Our main result is the following theorem.

Theorem 1**.**

Let us assume that $E$ is fixed. Let $\nu_{\min}$ and $\nu_{\max}$ be such that

[TABLE]

and

[TABLE]

Then, for all $\epsilon>0$, and for all $c_{O}>0$ and $c_{U}>c_{O}$ such that

[TABLE]

for any $R\in\{1,\ldots,n\}$ and for $c$ and $\rho$ sufficiently small universal constants, we have

[TABLE]

This theorem gives a prediction bound on the difference between the true and the estimated cluster matrices filtered by the matrix of means. Up to our knowledge, this is the first oracle bound for clustering using a generalized Bayesian NMF. Note that the oracle bound is not sharp as the leading constant is $1+\epsilon>1$, however $\epsilon$ may be chosen arbitrarily close to [math].

Note also that the claim that this result is PAC-Bayesian-flavored comes from the fact that the prediction machinery is largely inspired by Alquier and Guedj [2017], and the scheme of proof builds upon the PAC-Bayesian bound from Dalalyan and Tsybakov [2008]. Hence we kept that PAC-Bayesian filiation, even though the bound holds in expectation.

The following Theorem gives a more precise bound in the case where the noise $E$ is assumed to be iid Gaussian.

Theorem 2**.**

Assume that the dimension is larger that the number of clusters, i.e. $d>K$. Fix $\nu_{\max}>\nu_{\min}>0$ 111as the reader will be able to check, the values of $\nu_{\max}$ and $\nu_{\min}$ will play an essential role in the expression of the success probability of the method.. In addition to the assumptions of Theorem 1, assume that $E$ is iid Gaussian with minimum (resp. maximum) one-dimensional variance $\sigma_{\min}^{2}$ (resp. $\sigma_{\max}^{2}$) and assume also that the $\mu_{k}$ have Euclidean norm less that $1$ and pairwise scalar products less than $\mu$ in absolute value. Then, as long as $\mu<1/(K-1)$, for all $\epsilon>0$, and for all $c_{O}>0$ and $c_{U}>c_{O}$ such that

[TABLE]

with probability at least

[TABLE]

with

[TABLE]

we have

[TABLE]

with222Notice that $t_{\min}$ needs to be sufficiently smaller that $2/e$ in order for the term $K$ to become small for $n$ sufficiently large.

[TABLE]

and

[TABLE]

This theorem shows that a bound on the difference of the cluster matrices can be obtained when the matrix of means is sufficiently incoherent. Notice that this bound is not exactly component-wise, but considers a sum over clusters, which is perfectly relevant because the matrix $M$ does not distinguish between points in the same cluster. As expected, the smaller the coherence $\mu$, the better the oracle bound.

The proof of Theorem 2 is deferred to the supplementary material.

3 A Langevin sampler

In this section, we present a Langevin sampler for our estimator $\hat{U}_{\lambda},\hat{O}_{\lambda}$ 333Notation-wise, we will identify the Stiefel manifold with the set of matrices whose first $R$ columns form an orthonormal family and the remaining $n-R$ columns are set to zero.. Langevin-type samplers were first proposed by Grenander [1983], Grenander and Miller [1994], Roberts and Tweedie [1996], and have attracted a lot a attention lately in the statistical learning community Dalalyan [2017], Durmus and Moulines [2017], Brosse et al. [2018].

3.1 Computing the gradient on the Stiefel manifold

We start with some preliminary material about gradient computation on the Stiefel manifold from Edelman et al. [1998]. The Stiefel manifold can be interpreted as the set of equivalence classes

[TABLE]

As can easily be deduced from this quotient representation of the Stiefel manifold, the tangent space to the Stiefel manifold at a point $O$ is

[TABLE]

The canonical metric at a point $O$ is given by

[TABLE]

For $\Delta\in T_{O}(\mathbb{O}_{d,R})$, the exponential map is given by $O(t)=Oe^{t\Delta}I_{d,R}.$ The gradient at $O$ of a function $f$ defined on the Stiefel manifold $\mathbb{O}_{d,R}$ is given by 444This formula can be obtained using differentiation along the geodesic defined by the exponential map in the direction $\Delta$, for all $\Delta\in T_{O}(\mathbb{O}_{d,R})$.

[TABLE]

where $f_{O}(i,i^{\prime})=\frac{\partial f}{\partial\ O_{i,i^{\prime}}}$ for any $i$, $i^{\prime}=1,\ldots,n$.

3.2 The alternating Langevin sampler

The alternating Langevin sampler is described as follows. It consists in alternating between perturbed gradient steps in the matrix variable $O$ on the Stiefel manifold and perturbed gradient steps in the matrix variable $U$.

For clarity of the exposition, we give the formula for the gradient of $\|X-XUU^{t}\|_{F}^{2}$ as a function of $U$:

[TABLE]

Following Brosse et al. [2018], we propose the following algorithm.

In this algorithm the exponential function $\exp(O,H)$ at $O$ is given by [Edelman et al., 1998, Eq. 2.45] using different notation.

4 Numerical experiment

We provide some illustrations of the performance of our approach to clustering using a simple two component Gaussian Mixture and real satellite data. We performed an experiment with real high dimensional (satellite time series) data and one sample from the posterior after 300 iterations gave the result shown in Figure 1.

5 Proof of Theorem 2.1

We break down the proof in the following successive elementary steps.

5.1 Initial PAC-Bayesian bound

Theorem 3** (adapted from Dalalyan and Tsybakov, 2008).**

For $\lambda\leq 1/4$,

[TABLE]

where the infimum is taken over all probability measures $\rho$ which are absolutely continuous with respect to $\pi$. Here $KL$ denotes the Kullback-Leibler divergence.

5.2 Bounding the integral part

In order to bound the integral part in the bound given in (8), we define for any $R$ and any matrix $U^{0}\in\mathcal{U}_{R}\cap\mathbb{O}_{n,R,+}$ for any $c_{U},\ c_{O}\in(0,1]$, the measure

[TABLE]

Define $c=(c_{U},c_{O})$. Using these distributions we will be able to prove the following bound.

Lemma 4**.**

We have

[TABLE]

Proof.

Note that

[TABLE]

and thus,

[TABLE]

By the Cauchy-Schwartz inequality, we get

[TABLE]

Note further that since $(U,O)$ must belong to the support of $\rho_{R,U^{0},c_{U},c_{O}}$, we have

[TABLE]

Now, on the one hand, we have

[TABLE]

and since the columns of $U^{0}$ are orthonormal, we have $\|U^{0}\|=1$, which gives

[TABLE]

On the other hand,

[TABLE]

and using again that $\|U^{0}\|=1$,

[TABLE]

From this, we easily deduce that

[TABLE]

and

[TABLE]

which completes the proof. ∎

5.3 Upper bound on the Kullback-Leibler divergence

Lemma 5**.**

We have

[TABLE]

for some $\rho,c>0$ sufficiently small.

Proof.

By definition,

[TABLE]

We then have

[TABLE]

As long as $c_{U}^{2}\geq dR$, the inner probability can be bounded as follows [see equation 7.50 in Massart, 2007]:

[TABLE]

From this last inequality, we get

[TABLE]

We now use the elementary inequality $\log(1+x)\geq x/(1+x)$ for $x\in(-1,+\infty)$, and a lower bound on $\pi_{O}\left(\||O|-U^{0}\|_{F}^{2}\leq c_{O}^{2}\right)$ [given by Lee and Ruymgaart, 1996]. Therefore

[TABLE]

for $c$ and $\rho$ sufficiently small. From these two inequalities, we end the proof:

[TABLE]

∎

5.4 Combining the lemmæ

Combining the results of the previous lemmæ, we get the following proposition.

Proposition 6**.**

Let

[TABLE]

Then, for any $c_{U}$ such that $c_{U}\left(2+c_{U}\right)\leq\epsilon\frac{\nu}{\left(\|M\|+\|E\|\right)},$ we have

[TABLE]

Proof.

Let $\tilde{U}$ be a minimiser in the right hand-side in (10). Since

[TABLE]

and since $M=M\tilde{U}\tilde{U}^{t}$, we get $\|X-X\tilde{U}\tilde{U}^{t}\|_{F}=\|E(I-\tilde{U}\tilde{U}^{t})\|_{F}.$ As

[TABLE]

for all $U^{0}\in\mathbb{O}_{n,r,+}\cup\{U\mid MUU^{t}=M\}$, the claim then follows from combining the results of Lemma 4 and Lemma 5 above, and taking $c_{U}$ such that (3) in the main text holds. ∎

5.5 Assembling the elements

We have that

[TABLE]

Using Jensen’s inequality gives

[TABLE]

which, combined with Proposition 6 gives

[TABLE]

Given that

[TABLE]

this completes the proof of Theorem 2.1, since $U^{0}$ is any matrix satisfying the constraints.

6 Proof of Theorem 2.2

6.1 Control of $\mathbb{P}\left(\min_{U\in\mathbb{O}_{n,R}}\ \|E(I-UU^{t})\|_{F}\leq\nu_{\min}\right)$

The covering number of the Stiefel manifold in the operator norm was computed in Hinrichs et al. [2017] and is given by

[TABLE]

Let $\mathcal{N}_{\epsilon}$ denote an $\epsilon$-net in the operator norm for the Stiefel manifold with cardinality $N(\mathbb{O}_{n,R},\|\cdot\|,\epsilon)$. For any $U\in\mathbb{O}_{n,R}$, let $U^{\sharp}$ denote the closest matrix in $\mathcal{N}_{\epsilon}$ to $U$. Then, we have

[TABLE]

Moreover,

[TABLE]

Combining this last inequality with (12), we get

[TABLE]

Therefore, we obtain that

[TABLE]

where, for any $U$ in $\mathcal{N}_{\epsilon}$, $U^{\sharp}$ will denote the projection in operator norm of $U$ onto $\mathbb{O}_{n,R}$. Moreover

[TABLE]

Since $E$ is i.i.d. Gaussian with minimum one-dimensional variance $\sigma_{\min}^{2}$, for any $U^{\sharp}\in\mathbb{O}_{n,R}$, the lower tail of $\sigma_{\min}^{-1}\|E(I-U^{\sharp}U^{\sharp^{t}})\|_{F}^{2}$ is dominated by the lower tail of a $\chi^{2}(n(n-R))$ distribution and therefore, as recalled in Chrétien and Darses [2014, Lemma B1]

[TABLE]

Let us now tune $t$ and $\eta$. On the one hand, by Boucheron et al. [2013], we have that

[TABLE]

Thus, we will choose $\eta=\sigma\sqrt{dn+u}$. On the other hand, we will take $t$ such that

[TABLE]

i.e.

[TABLE]

Therefore, using the union bound we get

[TABLE]

6.2 Control of $\mathbb{P}\left(\max_{U\in\mathbb{O}_{n,R}}\ \|E(I-UU^{t})\|_{F}\geq\nu_{\max}\right)$

The same strategy applies, with slight modifications. We consider the same $\epsilon$-net as in the previous subsection. We first easily get

[TABLE]

Therefore,

[TABLE]

where, for any $U$ in $\mathcal{N}_{\epsilon}$, $U^{\sharp}$ again denotes the projection in operator norm of $U$ onto $\mathbb{O}_{n,R}$. We then get

[TABLE]

Since $E$ is i.i.d. Gaussian with maximum one-dimensional variance $\sigma_{\max}^{2}$, the upper tail of $\sigma_{\max}^{-1}\|E(I-U^{\sharp}U^{\sharp^{t}})\|_{F}^{2}$ is dominated by that of a $\chi^{2}(n(n-R))$ distribution for any $U^{\sharp}$ in $\mathbb{O}_{n,R}$, and therefore, as recalled in Chrétien and Darses [2014, Lemma B1]

[TABLE]

We can now tune $t$ and $\eta$. Recall that for all $u>0$, we have

[TABLE]

Thus, we will choose as before $\eta=\sigma_{\max}\sqrt{dn+u}$. On the other hand, we will take $t$ such that

[TABLE]

i.e.

[TABLE]

Therefore, using the union bound we get

[TABLE]

6.3 Control of $\|E\|$

We will also need to control $\|E\|$. Using Vershynin [2018, Section 4.4], we obtain

[TABLE]

with probability at least $1-\exp(-d)$.

6.4 Control of $\|M\|$

Finally, we need to compute $\|M\|$ as a function of the coherence $\mu$. Using the Gershgorin bound, we easily obtain

[TABLE]

6.5 End of the proof

Let us now study $\|M(T^{*}\hat{U}_{\lambda}\hat{U}_{\lambda}^{t})\|_{F}^{2}$. Let

[TABLE]

Let $\Upsilon^{\dagger}$ denote the pseudo inverse of $\Upsilon$, i.e.

[TABLE]

In particular, let $\pi$ denote a permutation which orders the data cluster wise, i.e. all data from cluster 1, followed by all data from cluster 2, …, all data from cluster $K$ and let $M_{\pi}$ denote the matrix whose columns are permuted by $\pi$. Then,

[TABLE]

Denote by $T^{*}_{\pi}$ the matrix obtained from $T^{*}$ after reordering its rows and columns using $\pi$, and $\hat{U}_{\lambda,\pi}$ denote the matrix obtained from $\hat{U}_{\lambda}$ after reordering its row using $\pi$.

[TABLE]

On the other hand,

[TABLE]

and we get that

[TABLE]

and the result follows.