# Optimal Clustering by Lloyd’s Algorithm for Low-Rank Mixture Model

Zhongyuan Lyu and Dong Xia111Dong Xia’s research was partially supported by
Hong Kong RGC Grant GRF 16300121 and GRF 16301622.
Department of Mathematics, Hong Kong University of Science and Technology

(())

###### Abstract

This paper investigates the computational and statistical limit in clustering
matrix-valued observations. We propose a low-rank mixture model (LrMM),
adapted from the classical Gaussian mixture model (GMM) to treat matrix-valued
observations, which assumes low-rankness for population center matrices. A
computationally efficient clustering method is designed by integrating Lloyd’s
algorithm and low-rank approximation. Once well-initialized, the algorithm
converges fast and achieves an exponential-type clustering error rate that is
minimax optimal. Meanwhile, we show that a tensor-based spectral method
delivers a good initial clustering. Comparable to GMM, the minimax optimal
clustering error rate is decided by the separation strength, i.e, the minimal
distance between population center matrices. By exploiting low-rankness, the
proposed algorithm is blessed with a weaker requirement on the separation
strength. Unlike GMM, however, the computational difficulty of LrMM is
characterized by the signal strength, i.e, the smallest non-zero singular
values of population center matrices. Evidences are provided showing that no
polynomial-time algorithm is consistent if the signal strength is not strong
enough, even though the separation strength is strong. Intriguing differences
between estimation and clustering under LrMM are discussed. The merits of low-
rank Lloyd’s algorithm are confirmed by comprehensive simulation experiments.
Finally, our method outperforms others in the literature on real-world
datasets.

##  1 Introduction

Nowadays, clustering matrix-valued observations becomes a ubiquitous task in
diverse fields. For instance, each highly variable region (HVR) in the var
genes of human malaria parasite (Larremore et al., 2013; Jing et al., 2021) is
representable by an adjacency matrix and a key scientific question is to
identify structurally-similar HVRs by, say, clustering the associated
adjacency matrices. The international trade flow of a commodity across
different countries can be viewed as a weighted adjacency matrix (Lyu et al.,
2021; Cai et al., 2022). Finding the similarity between the trading patterns
of different commodities is of great value in understanding the global
economic structure. This can also be achieved by clustering the weighted
adjacency matrices. Other notable examples include clustering multi-layer
social networks (Dong et al., 2012; Han et al., 2015) and multi-view data
(Kumar et al., 2011; Mai et al., 2021), modeling the connectivity of brain
networks (Arroyo et al., 2021; Sun and Li, 2019), clustering the correlation
networks between bacterial species (Stanley et al., 2016), and EEG data
analysis (Gao et al., 2021), etc.

Since matrix-valued observations can always be vectorized, a naive approach is
to ignore the matrix structure so that numerous classical clustering
algorithms, e.g. K-means or spectral clustering, are readily applicable.
However, matrix observations are usually blessed with hidden low-dimensional
structures, among which low-rankness is perhaps the most common and explored.
Network models such as stochastic block model (Holland et al., 1983; Jing et
al., 2021), random dot product graph (Athreya et al., 2017) and latent space
model (Hoff et al., 2002) often assume a low-rank expectation of adjacency
matrix. Low-rank structures have also been successfully explored in brain
image clustering (Sun and Li, 2019), EEG data analysis (Gao et al., 2021), and
international trade flow data (Lyu et al., 2021), to name but a few. Table 1
presents a summary of datasets analyzed in our paper, where the matrix ranks
$r_{k}$’s (suggested by the numerical performance of our algorithm) are much
smaller than the ambient dimensions $(d_{1},d_{2})$. Without loss of
generality, we assume $d_{1}\geq d_{2}$. For these applications, the naive
clustering approach becomes statistically sub-optimal since the planted low-
dimensional structure is overlooked.

Dataset | $n$ | $(d_{1},d_{2})$ | $K$ | Ranks
---|---|---|---|---
BHL (Mai et al., 2021) | 27 | (1124,4) | 3 | $\sim\\{1,1,1\\}$
EEG (Zhang et al., 1995) | 122 | (256,64) | 2 | $\sim\\{2,1\\}$
Malaria gene networks (Larremore et al., 2013) | 9 | (212,212) | 6 | $\leq 15$
UN trade flow networks (Lyu et al., 2021) | 97 | (48,48) | 2 | $\sim\\{3,2\\}$
Table 1: Summary of datasets. Here, $n$ is the sample size, $(d_{1},d_{2})$ is
the dimension of each matrix observation, and $K$ is number of clusters. The
underlying rank $(r_{k}^{\prime}s)$ of population center matrices from
different clusters can be unequal.

Motivated by the aforementioned applications, throughout this paper, we assume
that each matrix-valued observation has a low-rank expectation and the
expectations are equal for observations from the same cluster. It is the
essence of low-rank mixture model (LrMM), which shall be formally defined in
Section 2. Several clustering methods exploiting low-rankness have emerged in
the literature. Sun and Li (2019) introduces a tensor Gaussian mixture model
and recasts the clustering task as estimating the factors in low-rank tensor
decomposition. K-means clustering is then applied to the estimated factors.
While a sharp estimation error rate is derived under a suitable signal-to-
noise ratio (SNR) condition, the accuracy of clustering is not provided. A
tensor normal mixture model is proposed by Mai et al. (2021), where the
authors designed an enhanced EM algorithm for estimating the distributional
parameters. Under appropriate conditions, sharp estimation error rates are
established showing that minimax optimal test clustering error rate is
attainable. However, the training clustering error is missing, and it is even
unclear whether the proposed EM algorithm can consistently recover the true
cluster memberships. Aimed at analyzing multi-layer networks, Jing et al.
(2021) proposed a mixture multi-layer SBM where a spectral clustering method
based on tensor decomposition is investigated. Clustering error rate is
established under a fairly weak network sparsity condition, although the rate
is likely sub-optimal. More recently, Lyu et al. (2021) extended the mixture
framework to latent space model and a sub-optimal clustering error rate was
also derived. Note that Jing et al. (2021) and Lyu et al. (2021) both require
a rather restrictive condition in that $n=O(d_{1})$ rendering their theories
unattractive in many scenarios. Other representative works include Chen et al.
(2020), Cai et al. (2021), Gao et al. (2021) and Stanley et al. (2016), but
clustering error rates were not studied.

Note that LrMM reduces to the famous Gaussian mixture model (GMM) in the
dimension $d^{\ast}:=d_{1}d_{2}$ if each matrix-valued observation has a full-
rank expectation, and the noise matrix has i.i.d. standard normal entries.
Under GMM, Löffler et al. (2021) proved that a spectral method attains, with
high probability, an average mis-clustering error rate $\exp(-\Delta^{2}/8)$
that is optimal in the minimax sense. Here $\Delta$ is the minimal Euclidean
distance between the expected centers of distinct clusters (i.e., population
center matrices), referred to as the separation strength. This exponential
rate was established by Löffler et al. (2021) under a separation
strength222For narration simplicity, we set the number of clusters $K=O(1)$
here. condition $\Delta\gg 1+d^{\ast}n^{-1}$. Gao and Zhang (2022)
investigated a more general iterative algorithm that achieves the same
exponential rate under a weaker separation strength condition $\Delta\gg
1+(d^{\ast}/n)^{1/2}$. More recently, Zhang and Zhou (2022) applied the leave-
one-out method and proved the optimality of spectral clustering under a
relaxed separation strength condition. Besides deriving the optimal clustering
error rate, prior works also made efforts to establish the phase transitions
in exact recovery, i.e., when the clustering error is zero. Ndaoud (2018)
investigated a power iteration algorithm for a two-component GMM and proved
that exact recovery is attained w.h.p. if $\Delta^{2}$ is greater than
$\big{(}1+(1+2d^{\ast}n^{-1}\log^{-1}n)^{1/2}\big{)}\cdot\log n$. In addition,
the author showed that exact recovery is impossible if $\Delta^{2}$ is smaller
than the aforesaid threshold. Later, Chen and Yang (2021) established a
similar phase transition for general $K$-component GMM based on a semidefinite
programming (SDP) relaxation. These foregoing works suggest an intriguing gap
in the regime $n=O(d^{\ast})$: Ndaoud (2018) and Chen and Yang (2021) revealed
that exact recovery is achievable beyond the separation strength threshold
$(2d^{\ast}n^{-1}\log n)^{1/4}$, whereas the exponential-type clustering error
rate (Gao and Zhang, 2022; Zhang and Zhou, 2022) was derived only beyond the
threshold $(d^{\ast}/n)^{1/2}$. To our best knowledge, the gap still exists at
the moment. Jin et al. (2017) proposed a two-component symmetric sparse GMM
and investigated the phase transition in consistent clustering. Specifically,
they showed that, ignoring log factors, $\Delta\gg 1+s/n$ is necessary for
consistent clustering without restricting the computational complexity. Here
$s$ is the sparsity of the expected observation. A recent work (Löffler et
al., 2020) designed an SDP-based spectral method and established an
exponential-type clustering error rate when $\Delta$ is greater than
$1+s^{1/2}\log^{1/4}(d^{\ast})n^{-1/4}$. Moreover, they provided evidence
supporting the claim that no polynomial-time algorithm can consistent recover
the clusters if $\Delta$ is smaller than the aforesaid threshold, i.e., there
exists a statistical-to-computational gap for clustering in sparse GMM. Both
Jin et al. (2017) and Löffler et al. (2020) implied that the necessary
separation strength primarily depends on the intrinsic dimension $s$ rather
than the ambient dimension $d^{\ast}$. We remark that there is a vast
literature studying the clustering problem for GMM. A representative but
incomplete list includes Lu and Zhou (2016); Balakrishnan et al. (2017);
Dasgupta (2008); Fei and Chen (2018); Hajek et al. (2016); Verzelen and Arias-
Castro (2017); Witten and Tibshirani (2010); Abbe et al. (2020) and references
therein.

In contrast, the understanding of the limit of clustering for LrMM is still at
its infant stage. In this paper, we fill the void in the optimal clustering
error rate of LrMM and demonstrate that the rate is achievable by a
computationally fast algorithm. Challenges are posed from multiple fronts.
First of all, designing a computationally fast clustering procedure that
sufficiently exploits low-rank structure is non-trivial. Unlike (sparse) GMM
(Chen and Yang, 2021; Löffler et al., 2020), convex relaxation seems not
immediately accessible for the clustering of LrMM, especially when there are
more than two clusters. Non-convex approaches based on tensor decomposition
and spectral clustering (Jing et al., 2021; Luo and Zhang, 2022; Xia and Zhou,
2019) usually cannot distinguish the sample size dimension (i.e., $n$) and
data point dimension (i.e., $d_{1},d_{2}$). Their theoretical results become
sub-optimal when the sample size is much larger than $d_{1}$. On the technical
front, low-rankness makes deriving an exact exponential-type clustering error
rate even more difficult. Under GMM (Gao and Zhang, 2022; Löffler et al.,
2021), the exponential-type clustering error rate is established by carefully
studying the concentration phenomenon of a Gaussian linear form that usually
admits an explicit representation. Estimating procedures under LrMM, however,
often require multiple iterations of low-rank approximation, say, by singular
value decomposition (SVD). Consequently, deriving the concentration property
of respective linear forms under LrMM is much more involved than that under
GMM. Moreover, prior related works (Löffler et al., 2020; Jin et al., 2017;
Zhang and Xia, 2018; Lyu and Xia, 2022) provided evidences that imply the
existence of a statistical-to-computational gap. It is unclear which model
parameter characterizes such a gap and how the gap depends on the sample size
and dimensions. For instance, how the low-rankness benefits the separation
strength requirement? Interestingly, we discover that the gap is not
determined by the separation strength $\Delta$ but rather by the signal
strength (to be defined in Section 2) of the population center matrices.

Our main contributions are summarized as follows. First, we propose a
computationally fast clustering algorithm for LrMM. At its essence is the
combination of Lloyd’s algorithm (Lloyd, 1982; Lu and Zhou, 2016) and low-rank
approximation. Basically, given the updated cluster memberships of each
observation, the cluster centers are obtained by the SVD of the sample average
within each cluster. The whole algorithm involves only K-means clustering and
matrix SVDs. Secondly, we prove that, equipped with a good initial clustering,
the low-rank Lloyd’s algorithm converges fast and achieves the minimax optimal
clustering error rate $\exp(-\Delta^{2}/8)$ with high probability as long as
the separation strength satisfies $\Delta^{2}\gg 1+d_{1}r_{\textsf{\tiny
max}}/n$ and the signal strength is strong enough. Here $r_{\textsf{\tiny
max}}$ is the maximum rank among all the population center matrices. This
dictates that a weaker separation strength is sufficient for clustering under
LrMM if the rank $r_{\textsf{\tiny max}}=O(1)$. Our key technical tool to
develop the exponential-type error rate is a spectral representation formula
from Xia (2021), which has helped push forward the understanding of
statistical inference for low-rank models (Xia and Yuan, 2021; Xia et al.,
2022). Thirdly, we propose a novel tensor-based spectral method for obtaining
an initial clustering. Under similar separation strength and signal strength
conditions, this method delivers an initial clustering that is sufficiently
good for ensuring the convergence of low-rank Lloyd’s algorithm. Lastly,
compared with GMM that only requires a separation strength condition (Löffler
et al., 2020; Gao and Zhang, 2022), an additional signal strength condition
seems necessary under LrMM. We provide evidences, based on the low-degree
framework (Kunisky et al., 2019), showing that if the signal strength
condition fails, all polynomial-time algorithms cannot consistently recover
the true clusters, even when the separation strength is much stronger than the
aforesaid one. It is worth pointing out that, unlike tensor-based approaches
(Jing et al., 2021; Luo and Zhang, 2022; Xia and Zhou, 2019), our theoretical
results impose no constraints on the relation between $n$ and $(d_{1},d_{2})$.

The rest of the paper is organized as follows. Low-rank mixture model is
formalized in Section 2, and we introduce the low-rank Lloyd’s algorithm and a
tensor-based method for spectral initialization. The convergence performance
of Lloyds’ algorithm, minimax optimal exponential-type clustering error rate,
and guarantees of a tensor-based spectral initialization are established in
Section 3. We discuss the computational barriers of LrMM in Section 4. In
Section 5, we slightly modify the low-rank Lloyd’s algorithm and derive the
same minimax optimal clustering error rate requiring a slightly weaker signal
strength condition. We discuss the difference between estimation and
clustering under LrMM in Section 6. Further discussions are provided in
Section 7. Numerical simulations and real data examples are presented in
Section 8. All proofs and technical lemmas are relegated to the appendix.

##  2 Methodology

###  2.1 Background and notations

For nonnegative $D_{1},D_{2}$ , the notation $D_{1}\lesssim D_{2}$
(equivalently, $D_{2}\gtrsim D_{1}$) means that there exists an absolute
constant $C>0$ such that $D_{1}\leq CD_{2}$; $D_{1}\asymp D_{2}$ is equivalent
to $D_{1}\lesssim D_{2}$ and $D_{2}\lesssim D_{1}$, simultaneously. Let
$\|\cdot\|$ denote the $\ell_{2}$ norm for vectors and operator norm for
matrices, and $\|\cdot\|_{\rm F}$ denotes the matrix Frobenius norm. Denote
$\sigma_{1}(\mathbf{M})\geq\cdots\geq\sigma_{r}(\mathbf{M})>0$ the non-
increasing singular values of $\mathbf{M}$ where $r={\rm rank}(\mathbf{M})$.
We also define $\sigma_{\textsf{\tiny
min}}(\mathbf{M}):=\sigma_{r}(\mathbf{M})$. A third order tensor is a three-
dimensional array. Throughout the paper, a tensor is written in the
calligraphic bold font, e.g. ${\bm{\mathcal{M}}}\in\mathbb{R}^{d_{1}\times
d_{2}\times n}$. We use ${\mathscr{M}}_{1}({\bm{\mathcal{M}}})$ to denote the
mode-$1$ matricization of ${\bm{\mathcal{M}}}$ such that
${\mathscr{M}}_{1}({\bm{\mathcal{M}}})\in\mathbb{R}^{d_{1}\times(d_{2}n)}$ and
${\mathscr{M}}_{1}({\bm{\mathcal{M}}})(i_{1},(i_{2}-1)n+i_{3})={\bm{\mathcal{M}}}(i_{1},i_{2},i_{3}),\forall
i_{1}\in[d_{1}],i_{2}\in[d_{2}],i_{3}\in[n]$. The mode-$2$ and mode-$3$
matricizations are defined in a similar fashion. Then $\big{\\{}{\rm
rank}\big{(}{\mathscr{M}}_{k}({\bm{\mathcal{M}}})\big{)}:k=1,2,3\big{\\}}$ are
called Tucker rank or multilinear rank. The mode-$1$ marginal multiplication
between ${\bm{\mathcal{M}}}$ and a matrix
$\mathbf{U}^{\top}\in\mathbb{R}^{r\times d_{1}}$ results into a tensor of size
$r_{1}\times d_{2}\times n$, whose elements are

$\big{(}{\bm{\mathcal{M}}}\times_{1}\mathbf{U}^{\top}\big{)}(j_{1},i_{2},i_{3}):=\sum_{i_{1}=1}^{d_{1}}{\bm{\mathcal{M}}}(i_{1},i_{2},i_{3})\mathbf{U}(i_{1},j_{1}),\quad\forall
j_{1}\in[r],i_{2}\in[d_{2}],i_{3}\in[n]$

Similarly, we can define the mode-$2$ and mode-$3$ marginal multiplication.
Given ${\bm{\mathcal{S}}}\in\mathbb{R}^{r_{1}\times r_{2}\times
r_{3}},\mathbf{V}\in\mathbb{R}^{d_{2}\times
r_{2}},\mathbf{W}\in\mathbb{R}^{n\times r_{3}}$, the multi-linear product
${\bm{\mathcal{M}}}:={\bm{\mathcal{S}}}\times_{1}\mathbf{U}\times_{2}\mathbf{V}\times_{3}\mathbf{W}$
outputs a $d_{1}\times d_{2}\times n$ tensor defined by,

${\bm{\mathcal{M}}}(i_{1},i_{2},i_{3}):=\sum_{j_{1}=1}^{r_{1}}\sum_{j_{2}=1}^{r_{2}}\sum_{j_{3}=1}^{r_{3}}{\bm{\mathcal{S}}}(j_{1},j_{2},j_{3})\mathbf{U}(i_{1},j_{1})\mathbf{V}(i_{2},j_{2})\mathbf{W}(i_{3},j_{3})$
(1)

More details can be found in Kolda and Bader (2009). Denote $\mathbb{O}_{d,r}$
the set of all $d\times r$ matrices $\mathbf{U}$ such that
$\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}_{r}$, where $\mathbf{I}_{r}$ is the
$r\times r$ identity matrix. Eq. (1) is known as the Tucker decomposition if
$r_{k}={\rm rank}\big{(}{\mathscr{M}}_{k}({\bm{\mathcal{M}}})\big{)}$,
$\mathbf{U}\in\mathbb{O}_{d_{1},r_{1}},\mathbf{V}\in\mathbb{O}_{d_{2},r_{2}}$,
and $\mathbf{W}\in\mathbb{O}_{n,r_{3}}$.

###  2.2 Low-rank sub-Gaussian mixture

Suppose that the $d_{1}\times d_{2}$ matrix-valued observations
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ are i.i.d., and each of them has a
latent label $s_{i}^{\ast}\in[K]$. Here $K$ denotes the number of underlying
clusters, and without loss of generality, assume $d_{1}\geq d_{2}$. We assume
that there exists $K$ deterministic but unknown matrices
$\mathbf{M}_{1},\cdots,\mathbf{M}_{K}$ such that, conditioned on
$s_{i}^{\ast}=k$, $\mathbf{X}_{i}$ has i.i.d. zero-mean sub-Gaussian entries
with the mean matrix $\mathbf{M}_{k}$. This implies that
$\mathbf{X}_{i}|s_{i}^{\ast}=k$ is equal to $\mathbf{M}_{k}+\mathbf{E}_{i}$ in
distribution where the noise matrix $\mathbf{E}_{i}$ satisfies the following
assumption:

######  Assumption 1.

(Sub-Gaussian noise) The noise matrix $\mathbf{E}_{i}$ has i.i.d. zero-mean
entries and unit variance, and for
$\forall\mathbf{M}\in\mathbb{R}^{d_{1}\times d_{2}}$, the following
probability holds

$\mathbb{P}(\langle\mathbf{M},\mathbf{E}_{i}\rangle\geq t)\leq
e^{-t^{2}/(2\sigma_{\textsf{sg}}^{2}\cdot\|\mathbf{M}\|_{\rm
F}^{2})},\quad\forall t>0,$

where $\sigma_{\textsf{sg}}>0$ is the sub-Gaussian constant.

Throughout the paper, we let $\sigma_{\textsf{sg}}^{2}=1$ without loss
generality (say, by substituting $\mathbf{X}_{i}$ with
$\mathbf{X}_{i}/\sigma_{\textsf{sg}}$). Moreover, we assume that the latent
labels $s_{1}^{\ast},\cdots,s_{n}^{\ast}$ are i.i.d. and

$\mathbb{P}(s_{i}^{\ast}=k)=\pi_{k},\quad\forall k\in[K];\quad{\rm
where}\quad\sum_{k=1}^{K}\pi_{k}=1.$ (2)

Here the unknown $\pi_{k}>0$ stands for the mass of $k$-th cluster. Put it
differently, the matrix-valued observations have a marginal distribution

$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\ \stackrel{{\scriptstyle{\rm
i.i.d.}}}{{\sim}}\ \sum_{k=1}^{k}\pi_{k}\cdot
p_{\mathbf{M}_{k},\sigma_{\textsf{sg}}^{2}}(\mathbf{X})$ (3)

where $p_{\mathbf{M}_{k},\sigma_{\textsf{sg}}^{2}}(\mathbf{X})$ is the density
function of matrix observation $\mathbf{X}\in\mathbb{R}^{d_{1}\times d_{2}}$
with independent entries of unit variance, the sub-Gaussian constant
$\sigma_{\textsf{sg}}$ and mean matrix $\mathbf{M}$. Let $r_{k}={\rm
rank}(\mathbf{M}_{k})$ and assume $r_{k}\ll d_{2}$ for all $k$, i.e., all the
population center matrices are low-rank. Model (3) is referred to as the low-
rank mixture model (LrMM). For simplicity, we treat the ranks $r_{k}$’s as
known and will briefly discuss how to estimate them in Section 7. We denote
the compact SVD of population center matrices by
$\mathbf{M}_{k}=\mathbf{U}_{k}\mathbf{\Sigma}_{k}\mathbf{V}_{k}^{\top}$ with
$\mathbf{U}_{k}\in\mathbb{O}_{d_{1},r_{k}}$ and
$\mathbf{V}_{k}\in\mathbb{O}_{d_{2},r_{k}}$. The signal strength of
$\mathbf{M}_{k}$ is characterized by $\sigma_{\textsf{\tiny
min}}(\mathbf{M}_{k}):=\sigma_{r_{k}}(\mathbf{M}_{k})$. We remark that
estimating $K$ is a challenging question even under GMM. Hence, throughout
this paper, it is assumed that $K$ is provided beforehand.

Sun and Li (2019) introduced a tensor Gaussian mixture model without
specifically imposing low-rank structures on the center matrices. A similar
tensor normal mixture model without low-rank assumptions is proposed by Mai et
al. (2021). Our LrMM can be viewed as a generalization of mixture multi-layer
SBM proposed by Jing et al. (2021) and as an extension of the symmetric two-
component case introduced by Lyu et al. (2021). Mixture of low-rank matrix
normal models have also appeared in Gao et al. (2021) for image analysis.

Since our goal of current paper is to investigate the fundamental limits of
clustering matrix-valued observations, hereafter, we view the latent labels
$s_{i}^{\ast},i\in[n]$ as a fixed realization sampled from the mixture
distribution (2). Then the matrix-valued observations can be written in the
following form:

$\mathbf{X}_{i}=\mathbf{M}_{s_{i}^{\ast}}+\mathbf{E}_{i},\quad i\in[n]$ (4)

Denote $\mathbf{s}^{\ast}=(s^{\ast}_{1},\cdots,s^{\ast}_{n})$ the collection
of true latent labels, known as the cluster membership vector. The size of
each cluster is given by
$n_{k}^{\ast}:=\sum_{i=1}^{n}\mathbb{I}\left(s_{i}^{\ast}=k\right),\forall
k\in[K].$ With mild conditions under LrMM, Chernoff bound (Chernoff, 1952)
guarantees $n_{k}^{\ast}\asymp n\pi_{k}$ with high probability.

Given an estimated cluster membership vector
$\widehat{}\mathbf{s}:=(\widehat{s}_{1},\cdots,\widehat{s}_{n})\in[K]^{n}$,
its clustering error is measured by the Hamming distance defined by

$h_{\textsf{c}}(\widehat{}\mathbf{s},\mathbf{s}^{\ast})=\min_{\pi:\text{
permutation of }[K]}\
\sum_{i=1}^{n}\mathbb{I}\left(\widehat{s}_{i}\neq\pi(s_{i}^{\ast})\right)$ (5)

For technical convenience, we also define the the following Frobenious error
related to ${\bm{\mathcal{M}}}$:

$\ell_{\textsf{c}}(\widehat{}\mathbf{s},\mathbf{s}^{\ast})=\min_{\pi:\text{
permutation of }[K]}\
\sum_{i=1}^{n}\big{\|}\mathbf{M}_{\widehat{s}_{i}}-\mathbf{M}_{\pi(s_{i}^{\ast})}\big{\|}_{\rm
F}^{2}.$

###  2.3 Low-rank Lloyd’s algorithm

Lloyd’s algorithm (Lloyd, 1982) or K-means algorithm is perhaps, conceptually
and implementation-wise, the most simple yet effective method for clustering.
It is an iterative algorithm, which consists of two main routines at each
iteration: 1). provided with an estimated cluster membership vector, the
cluster centers are updated by taking the sample average within every
estimated cluster; 2). provided with the updated cluster centers, every data
point is assigned an updated cluster label according to its distances from the
cluster centers. The iterations are terminated once converged. The success of
Lloyd’s algorithm is highly reliant on a good initial clustering or initial
cluster centers. It is proved by Lu and Zhou (2016) and Gao and Zhang (2022)
that, if well initialized, Lloyd’s algorithm converges fast and achieves
minimax optimal clustering error for GMM and community detections under
stochastic block model.

The original Lloyd’s algorithm updates the cluster centers by taking the
vanilla sample average. This approach is sub-optimal under LrMM because the
underlying low-rank structure is overlooked. It is well-known that exploiting
the low-rankness can further de-noise the estimates. Towards that end, we
propose the low-rank Lloyd’s algorithm whose details are enumerated in
Algorithm 1. Compared with the original Lloyd’s algorithm, the low-rank
version only modifies the procedure of updating the cluster centers. At the
$(t+1)$-th iteration, given the current cluster labels
$\widehat{}\mathbf{s}^{(t)}$ and for each $k$, we calculate the sample average
$\bar{}\mathbf{X}_{k}(\widehat{}\mathbf{s}^{(t)})$ defined as in Algorithm 1,
and then update the cluster center by

$\widehat{}\mathbf{M}_{k}^{(t+1)}:=\widehat{}\mathbf{U}_{k}^{(t)}\widehat{}\mathbf{U}_{k}^{(t)\top}\bar{}\mathbf{X}_{k}(\widehat{}\mathbf{s}^{(t)})\widehat{}\mathbf{V}_{k}^{(t)}\widehat{}\mathbf{V}_{k}^{(t)\top}$

where $\widehat{}\mathbf{U}_{k}^{(t)}$ and $\widehat{}\mathbf{V}_{k}^{(t)}$
are the top-$r_{k}$ left and right singular vectors of
$\bar{\mathbf{X}}_{k}(\widehat{}\mathbf{s}^{(t)})$, respectively. The update
of cluster labels is unchanged compared with the original Lloyd’s algorithm.

Algorithm 1 Low-rank Lloyd’s Algorithm (lr-Lloyd)

Input: Observations
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\in\mathbb{R}^{d_{1}\times d_{2}}$,
initial estimate $\widehat{}\mathbf{s}^{(0)}$, ranks $\\{r_{k}\\}_{k=1}^{K}$.

for $t=1,\ldots,T$ do

for each $k=1,\cdots,K$: (update cluster centers)
$\widehat{}\mathbf{M}_{k}^{(t)}\leftarrow\textrm{best rank-}r_{k}\textrm{
approximation of
}\bar{\mathbf{X}}_{k}(\widehat{}\mathbf{s}^{(t-1)}):=\frac{\sum_{i=1}^{n}\mathbb{I}\left(\widehat{s}_{i}^{(t-1)}=k\right)\mathbf{X}_{i}}{\sum_{i=1}^{n}\mathbb{I}\left(\widehat{s}_{i}^{(t-1)}=k\right)}$
(6)

for each $i=1,\cdots,n$: (update cluster labels)
$\widehat{s}_{i}^{(t)}\leftarrow\underset{k\in[K]}{\arg\min}\
\|\mathbf{X}_{i}-\widehat{}\mathbf{M}_{k}^{(t)}\|_{\rm F}^{2}$

end for

Output: $\widehat{}\mathbf{s}:=\widehat{}\mathbf{s}^{(T)}$

Conceptually, our low-rank Lloyd’s algorithm is a direct adaptation of Lloyd’s
algorithm to accommodate low-rankness. However, the low-rank update of cluster
centers poses fresh and highly non-trivial challenges in studying the
convergence behavior of Algorithm 1. The original Lloyd’s algorithm simply
takes the sample average and thus admits a clean and explicit representation
form for the updated centers, which plays a critical role in technical
analysis, as in Gao and Zhang (2022). In sharp contrast, the required SVD in
Algorithm 1 involves intricate and non-linear operations on the matrix-valued
observations, and there is surely no clean and explicit representation form
for $\widehat{}\mathbf{M}_{k}^{(t)}$. More advanced tools are in need for our
purpose, as shall be explained in Section 3.

###  2.4 Tensor-based spectral initialization

The success of Algorithm 1 crucially depends on a reliable initial clustering.
A naive approach is to vectorize the matrix observations, concatenate them
into a new matrix of size $n\times(d_{1}d_{2})$, then borrow the classic
spectral clustering method as in Löffler et al. (2021) and Zhang and Zhou
(2022). Unfortunately, the naive approach turns out to be sub-optimal for
ignoring the planted low-dimensional structure in the row space.

Our proposed initial clustering is based on tensor decomposition. Towards that
end, we construct a third-order data tensor
${\bm{\mathcal{X}}}\in\mathbb{R}^{d_{1}\times d_{2}\times n}$ by stacking the
matrix-valued observations slice by slice, i.e., its $i$-th slice333We follow
Matlab syntax tradition and denote ${\cal X}(:,:,i)$ the sub-tensor by fixing
one index. ${\bm{\mathcal{X}}}(:,:,i)=\mathbf{X}_{i}$. The noise tensor
$\bm{\mathcal{E}}$ is defined in the same fashion. The signal tensor
${\bm{\mathcal{M}}}$ is constructed such that
${\bm{\mathcal{M}}}(:,:,i)=\mathbf{M}_{s_{i}^{\ast}}$. The tensor form of LrMM
(4) is

${\bm{\mathcal{X}}}={\bm{\mathcal{M}}}+\bm{\mathcal{E}}$ (7)

Interestingly, eq. (7) coincides with the famous tensor SVD or PCA model
(Zhang and Xia, 2018; Xia and Zhou, 2019; Liu et al., 2022). Let
$\mathring{r}:=\sum_{k=1}^{K}r_{k}$. Indeed, the signal tensor
${\bm{\mathcal{M}}}$ admits the following low-rank decomposition

${\bm{\mathcal{M}}}={\bm{\mathcal{S}}}\times_{1}\mathbf{U}\times_{2}\mathbf{V}\times_{3}\mathbf{W}$
(8)

where the $\mathring{r}\times\mathring{r}\times K$ core tensor
${\bm{\mathcal{S}}}$ is constructed as

${\bm{\mathcal{S}}}(:,:,k):={\rm diag}({\bf 0}_{r_{1}},\cdots,{\bf
0}_{r_{k-1}},\mathbf{\Sigma}_{r_{k}},{\bf 0}_{r_{k+1}},\cdots,{\bf
0}_{r_{K}})$

and
$\mathbf{U}=(\mathbf{U}_{1},\cdots,\mathbf{U}_{K})\in\mathbb{R}^{d_{1}\times\mathring{r}}$,
$\mathbf{V}=(\mathbf{V}_{1},\cdots,\mathbf{V}_{K})\in\mathbb{R}^{d_{2}\times\mathring{r}}$,
$\mathbf{W}=(\mathbf{e}_{s_{1}^{\ast}},\cdots,\mathbf{e}_{s_{n}^{\ast}})^{\top}\in\\{0,1\\}^{n\times
K}$. Here $\mathbf{e}_{k}$ denotes the $k$-th canonical basis vector in
Euclidean space whose dimension might vary at different appearances. Clearly,
the rows of $\mathbf{W}$ provide the cluster information and is referred to as
the cluster membership matrix. Note that (8) is not necessarily the Tucker
decomposition since $\mathbf{U},\mathbf{V}$ might be rank-deficient, in which
case the decomposition in the form (8) is not unique and
$\mathbf{U},\mathbf{V}$ become unrecoverable.

The singular space of ${\bm{\mathcal{M}}}$ is uniquely characterized by its
Tucker decomposition. To this end, denote
$\mathbf{U}^{\ast}\in\mathbb{O}_{d_{1},r_{\mathbf{U}}}$ and
$\mathbf{V}^{\ast}\in\mathbb{O}_{d_{2},r_{\mathbf{V}}}$ the left singular
vectors of $\mathbf{U}$ and $\mathbf{V}$, respectively. Here, $r_{\mathbf{U}}$
and $r_{\mathbf{V}}$ are the ranks of ${\mathscr{M}}_{1}({\bm{\mathcal{M}}})$
and ${\mathscr{M}}_{2}({\bm{\mathcal{M}}})$, respectively. Define
$\mathbf{W}^{\ast}\in\mathbb{O}_{n,K}$ by normalizing the columns of
$\mathbf{W}$. Re-compute the core tensor
${\bm{\mathcal{S}}}^{\ast}:={\bm{\mathcal{M}}}\times_{1}\mathbf{U}^{\ast\top}\times_{2}\mathbf{V}^{\ast\top}\times_{3}\mathbf{W}^{\ast\top}$
that is of size $r_{\mathbf{U}}\times r_{\mathbf{V}}\times K$. Finally, we re-
parameterize the signal tensor via its Tucker decomposition

${\bm{\mathcal{M}}}={\bm{\mathcal{S}}}^{\ast}\times_{1}\mathbf{U}^{\ast}\times_{2}\mathbf{V}^{\ast}\times_{3}\mathbf{W}^{\ast}$
(9)

Here $\mathbf{U}^{\ast},\mathbf{V}^{\ast},\mathbf{W}^{\ast}$ are usually
called the singular vectors of ${\bm{\mathcal{M}}}$. Still, the rows of
$\mathbf{W}^{\ast}$ tell the cluster information in that
$\mathbf{W}^{\ast}(i,:)=\mathbf{W}^{\ast}(j,:)$ iff
$s_{i}^{\ast}=s_{j}^{\ast}$, i.e, $i,j$ belongs to the same cluster. We note
that there are interesting special cases concerning the values of
$r_{\mathbf{U}},r_{\mathbf{V}}$. For instance, if
$r_{\mathbf{U}}=r_{\mathbf{V}}=r_{1}$, it implies that all the population
center matrices share the same low-dimensional singular space with
$\mathbf{M}_{1}$, which simplifies theoretical investigate of our proposed
initialization method. Another special case is
$r_{\mathbf{U}}=r_{\mathbf{V}}=\mathring{r}$, namely the singular spaces of
all population center matrices are separated to a certain degree. Intuitively,
the clustering problem becomes easier. See Section 3.2 for discussions of both
cases.

We now present our tensor-based spectral method for initial clustering. Unlike
the aforementioned naive spectral method, ours is specifically designed to
exploit the low-rank structure of ${\bm{\mathcal{M}}}$ in the $1$st and $2$nd
dimension. Without loss of generality, we treat $r_{\mathbf{U}}$ and
$r_{\mathbf{V}}$ as known here and shall discuss ways to estimate them in
Section 7. Our method consists of three crucial steps with details in
Algorithm 2. Step 1 aims to estimate the singular vectors $\mathbf{U}^{\ast}$
and $\mathbf{V}^{\ast}$. Here, higher order SVD (HOSVD) is obtained by
applying SVD to the matricizations ${\mathscr{M}}_{1}({\bm{\mathcal{M}}})$ and
${\mathscr{M}}_{2}({\bm{\mathcal{M}}})$. See, for instance, De Lathauwer et
al. (2000) and Xia and Zhou (2019). The estimated singular vectors are used
for denoising in Step 2 by projecting the noise into a low-dimensional space.
Step 3 applies the classical K-means clustering (Löffler et al., 2021; Zhang
and Zhou, 2022) to the denoised observations. Note that solving K-means is
generally NP-hard (Mahajan et al., 2009), but there exist fast algorithms
(Kumar et al., 2004) achieving an approximate solution.

Algorithm 2 Tensor-based Spectral Initialization (TS-Init)

Input: Observations
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\in\mathbb{R}^{d_{1}\times d_{2}}$ or a
tensor ${\bm{\mathcal{X}}}\in\mathbb{R}^{d_{1}\times d_{2}\times n}$ by
concatenating the matrix observations slice by slice.

  1. 1.

Obtain the estimated factor matrices $\widehat{}\mathbf{U}$ and
$\widehat{}\mathbf{V}$ by applying HOSVD to the tensor ${\bm{\mathcal{X}}}$ in
mode-$1$ and mode-$2$ with rank $r_{\mathbf{U}}$ and $r_{\mathbf{V}}$,
respectively.

  2. 2.

Project ${\bm{\mathcal{X}}}$ onto the column space of $\widehat{}\mathbf{U}$
and $\widehat{}\mathbf{V}$ by

$\widehat{\bm{\mathcal{G}}}:={\bm{\mathcal{X}}}\times_{1}\widehat{}\mathbf{U}\widehat{}\mathbf{U}^{\top}\times_{2}\widehat{}\mathbf{V}\widehat{}\mathbf{V}^{\top}\in\mathbb{R}^{d_{1}\times
d_{2}\times n}$

  3. 3.

Apply k-means on rows of
$\widehat{}\mathbf{G}:={\mathscr{M}}_{3}(\widehat{\bm{\mathcal{G}}})\in\mathbb{R}^{n\times
d_{1}d_{2}}$ to obtain initializer for $\mathbf{s}^{\ast}$, i.e.

$(\widehat{}\mathbf{s}^{(0)},\\{\widehat{}\mathbf{M}_{k}^{(0)}\\}_{k=1}^{K}):=\operatorname*{arg\,min}_{\mathbf{s}\in[K]^{n},\\{\mathbf{M}_{k}\\}_{k=1}^{K},\mathbf{M}_{k}\in\mathbb{R}^{d_{1}\times
d_{2}},\forall
k}\sum_{i=1}^{n}\left\|[\widehat{}\mathbf{G}]_{i\cdot}-vec(\mathbf{M}_{s_{i}})\right\|^{2}$

Output: $\widehat{}\mathbf{s}^{(0)}$

Algorithm 2 improves the naive spectral clustering whenever
$\widehat{}\mathbf{U}$ and $\widehat{}\mathbf{V}$ are reliable estimates of
their population counterparts. This suggests that a certain signal strength
condition on ${\mathscr{M}}_{1}({\bm{\mathcal{S}}}^{\ast})$ and
${\mathscr{M}}_{2}({\bm{\mathcal{S}}}^{\ast})$ is necessary. We remark that
the higher order orthogonal iteration (HOOI, Zhang and Xia (2018)) algorithm
for tensor decomposition is not suitable for our purpose since it requires a
lower bound on $\sigma_{\textsf{\tiny
min}}\big{(}{\mathscr{M}}_{3}({\bm{\mathcal{S}}}^{\ast})\big{)}$, which is too
restrictive under LrMM. See Section 3.2 for more explanations.

##  3 Minimax Optimal Clustering Error Rate of LrMM

In this section, we establish the convergence performance of low-rank Lloyd’s
algorithm, validate our tensor-based spectral initialization, and derive the
minimax optimal clustering error rate for LrMM (3). The hardness of clustering
under LrMM is determined primary by two quantities:

Separation strength $\displaystyle\Delta:=\min_{a\neq
b,a,b\in[K]}\left\|\mathbf{M}_{a}-\mathbf{M}_{b}\right\|_{\rm F}$

The separation strength is a generalization of the minimum $\ell_{2}$ distance
between different population centers under GMM (Lu and Zhou, 2016; Chen and
Yang, 2021; Gao and Zhang, 2022), which characterizes the intrinsic difficult
in clustering the observations. In fact, the minimax optimal error rate, i.e,
the best achievable clustering accuracy, is exclusively decided by $\Delta$.

###  3.1 Iterative convergence of low-rank Lloyd’s algorithm

The performance of Lloyd’s algorithm also relies on the minimal cluster size
(Lu and Zhou, 2016). To this end, define
$\alpha:=\min_{k\in[K]}n_{k}^{\ast}\cdot(n/K)^{-1}$, where recall that
$n_{k}^{\ast}:=|\\{i\in[n]:s_{i}^{\ast}=k\\}|$ is the size of $k$-th cluster.
The cluster sizes are said to be balanced if $\alpha\asymp 1$. The hamming
distance $h_{\textsf{c}}(\widehat{}\mathbf{s},\mathbf{s}^{\ast})$ is defined
as in eq. (5). Without loss of generality, we assume $r:=r_{1}$ is the largest
amongst $\\{r_{k}:k\in[K]\\}$ and $d:=d_{1}\geq d_{2}$.

Due to technical reasons, we define
$\kappa_{0}:={\max_{k\in[K]}\|\mathbf{M}_{k}\|}/\min_{k\in[K]}\sigma_{\textsf{\tiny
min}}(\mathbf{M}_{k})$, which can be viewed as the maximum condition number of
all population center matrices. It usually does not appear in the literature
of GMM, but is of unique importance under LrMM. This quantity plays a critical
role in connecting the accuracy of updated center matrix
$\widehat{}\mathbf{M}_{k}^{(t)}$ to the current clustering accuracy
$h_{\textsf{c}}(\widehat{}\mathbf{s}^{(t-1)},\mathbf{s}^{\ast})$. Since
$\widehat{}\mathbf{M}_{k}^{(t)}$ stems from the SVD of
$\bar{}\mathbf{X}_{k}(\widehat{}\mathbf{s}^{(t-1)})$, whose accuracy is
characterized by the strength of signal $\mathbf{M}_{k}$ and size of
perturbation
$\bar{}\mathbf{X}_{k}(\widehat{}\mathbf{s}^{(t-1)})-\mathbf{M}_{k}$. Besides
random noise, the latter term, roughly, consists of
$(n_{a}^{\ast})^{-1}h_{\textsf{c}}(\widehat{}\mathbf{s}^{(t-1)},\mathbf{s}^{\ast})\big{(}\mathbf{M}_{k^{\prime}\neq
k}-\mathbf{M}_{k}\big{)}$, whose operator norm can be controlled by
$O\big{(}(n_{a}^{\ast})^{-1}h_{\textsf{c}}(\widehat{}\mathbf{s}^{(t-1)},\mathbf{s}^{\ast})\kappa_{0}\sigma_{\textsf{\tiny
min}}(\mathbf{M}_{k})\big{)}$. Hence $\kappa_{0}$ is, perhaps, the unavoidable
price to be paid for taking advantage of low-rankness (chen2021learning).

The following theorem presents the convergence performance of low-rank Lloyd’s
algorithm (Algorithm 1). Due to the local nature of Lloyd’s algorithm, its
success highly relies on a good initialization. Theorem 1 assumes the initial
clustering is consistent, i.e., initial clustering error approaches zero
asymptotically as $n\to\infty$. Under suitable conditions of separation
strength and signal strength, the output of Algorithm 1 attains an
exponential-type error rate. The constant factor $1/8$ in the exponential rate
exactly matches the minimax lower bound in Theorem 3. Notice that our result
is non-asymptotic, and all asymptotic conditions in Theorem 1 are to guarantee
the sharp constant $1/8$ in eq. (12). More precisely, through a careful
inspection on our analysis, the implicit term $o(1)$ in the exponential rate
in Theorem 1 can be chosen at the order $\Omega\left(\left(Kr(d+\log n)(\alpha
n)^{-1}/\Delta^{2}\right)^{1/2-\epsilon}\right)=o(1)$ for any fixed
$\epsilon\in(0,1/2)$.

######  Theorem 1.

Suppose $d\geq C_{0}\log K$ for some absolute constant $C_{0}>0$. Assume that

  1. (i)

initial clustering error:

$n^{-1}\cdot\ell_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})=o\left(\frac{\alpha}{\kappa_{0}^{2}K}\Delta^{2}\right)$
(10)

  2. (ii)

separation strength:

$\frac{\Delta^{2}}{\alpha^{-1}(\kappa_{0}^{2}\vee
Kr)Kr\left(\frac{d}{n}+1\right)}\rightarrow\infty$ (11)

Let $\widehat{}\mathbf{s}^{(t)}$ be the cluster labels at $t$-th iteration
generated by Algorithm 1. Then, for all $t\geq 1$, we have

$n^{-1}\cdot
h_{\textsf{c}}(\widehat{}\mathbf{s}^{(t)},\mathbf{s})\leq\exp\left(-(1-o(1))\frac{\Delta^{2}}{8}\right)+\frac{1}{2^{t}}$
(12)

with probability at least $1-\exp(-\Delta)-\exp(-c_{0}d)$ with some absolute
constant $c_{0}>0$.

By Theorem 1, after at most $O\big{(}\min\\{\Delta^{2},\log n\\}\big{)}$
iterations, our low-rank Lloyd’s algorithm achieves the minimax optimal
clustering error rate $\exp(-\Delta^{2}/8)$, which is the same optimal rate
for classical GMM (Lu and Zhou, 2016; Löffler et al., 2020; Gao and Zhang,
2022; Zhang and Zhou, 2022) and is exclusively decided by the separation
strength $\Delta$. It is worth noting that provided with good initialization,
lr-Lloyd solely requires separation strength strong enough to achieve such
optimal rate.

Blessing of low-rankness and comparison with GMM. If low-rankness is ignored
so that LrMM is treated as GMM, the exponential-type error rate is established
only in the regime of separation strength $\Delta\gg 1+(d_{1}d_{2}/n)^{1/2}$
(Gao and Zhang, 2022; Zhang and Zhou, 2022). In contrast, our condition (11)
only requires $\Delta\gg 1+(d_{1}/n)^{1/2}$ if $r,K,\kappa_{0},\alpha=O(1)$.

Discussions on separation strength $\Delta$. The separation strength condition
is typical in the literature of clustering problems (Vempala and Wang, 2004;
Löffler et al., 2021). To see why our condition (11) is minimal, without loss
of generality, consider the case $\alpha\asymp 1$ and $K=2$. Moreover, assume
the singular vectors $\mathbf{U}_{1}=\mathbf{U}_{2}$ and
$\mathbf{V}_{1}=\mathbf{V}_{2}$, and they are already known. One can multiply
each observation by $\mathbf{U}_{1}^{\top}$ from left and by $\mathbf{V}_{1}$
from right, which reduces LrMM to GMM in the dimension $r^{2}$. Literature of
GMM (Gao and Zhang, 2022; Löffler et al., 2021; Zhang and Zhou, 2022) all
impose a separation strength condition $\Delta\gg 1$. This certifies the
constant $1$ in eq. (11). To understand the term $(rd/n)^{1/2}$, consider that
the true labels of first $n-1$ observations are revealed to us and our goal is
to estimate the label of the $n$-th sample $\mathbf{X}_{n}$. A natural way is
to first estimate the population centers utilizing the given labels
$\mathbf{s}^{\ast}_{1},\cdots,\mathbf{s}_{n-1}^{\ast}$, denoted by
$\widehat{}\mathbf{M}_{1}$ and $\widehat{}\mathbf{M}_{2}$, respectively. The
literature of matrix denoising (Cai and Zhang, 2018; Xia, 2021; Gavish and
Donoho, 2017) tells that the minimax optimal estimation error is at the order
$\|\widehat{}\mathbf{M}_{1}-\mathbf{M}_{1}\|_{\rm
F}\asymp\|\widehat{}\mathbf{M}_{2}-\mathbf{M}_{2}\|\asymp(rd/n)^{1/2}$. Thus
$\Delta\gg(rd/n)^{1/2}$ is necessary for consistently distinguishing the two
clusters. The above rationale suggests that our separation strength condition
(11) might be minimal up to the order of $n$, if only the exponential-type
error rate is sought.

We explained a gap concerning the separation strength in existing literature
of GMM. Under GMM with dimension $d^{\ast}=d_{1}d_{2}$ and $n\leq d^{\ast}$,
the exponential-type rate (Gao and Zhang, 2022; Zhang and Zhou, 2022) is
established in the regime $\Delta\gg(d^{\ast}/n)^{1/2}$, whereas exact
clustering results (Ndaoud, 2018; Chen and Yang, 2021) are attained in the
regime $\Delta\gtrsim(d^{\ast}n^{-1}\log n)^{1/4}$. This leaves a natural
question under LrMM: is the separation strength condition (11) is relaxable to
the scale $n^{-1/4}$? Unfortunately, answering this question is perhaps more
challenging than that under GMM. We note that Ndaoud (2018) and Chen and Yang
(2021) achieve the $O(n^{-1/4})$ barrier by focusing entirely on clustering
and by circumventing the estimation of population centers. Nonetheless, under
LrMM, exploiting the low-rank structure demands estimating the population
center matrices. We suspect, together with the aforementioned special
examples, that condition (11) might not be improvable in terms of the order of
$n$. Anyhow, It’s unclear whether one can obtain a sharper characterization of
$\Delta$ under LrMM using other methods like SDP. Further investigation in
this respect is out of the scope of current paper.

###  3.2 Guaranteed initialization

Besides the separation strength condition, Theorem 1 requires a consistent
initial clustering. We now demonstrate the validity of tensor-based Algorithm
2. Observe that denoising by spectral projection (Step 2 of Algorithm 2) is
only beneficial if $\widehat{}\mathbf{U}$ and $\widehat{}\mathbf{V}$ are
properly aligned with $\mathbf{U}^{\ast}$ and $\mathbf{V}^{\ast}$,
respectively. For that purpose, the signal strengths of
${\mathscr{M}}_{1}({\bm{\mathcal{M}}})$ and
${\mathscr{M}}_{2}({\bm{\mathcal{M}}})$, i.e., $\sigma_{\textsf{\tiny
min}}\big{(}{\mathscr{M}}_{1}({\bm{\mathcal{M}}})\big{)}$ and
$\sigma_{\textsf{\tiny
min}}\big{(}{\mathscr{M}}_{2}({\bm{\mathcal{M}}})\big{)}$, needs to be
sufficiently strong. For simplicity, we let $\Lambda_{\textsf{\tiny
min}}:=\min_{j=1,2}\left\\{\sigma_{\min}({\mathscr{M}}_{j}({\bm{\mathcal{M}}}))\right\\}$
denote tensor signal strength in 1st and 2nd modes of ${\bm{\mathcal{M}}}$, or
simply the tensor signal strength of ${\bm{\mathcal{M}}}$. Note that this is a
slightly different definition from classical tensor literature, where the
signal strength is usually defined as
$\min_{j=1,2,3}\left\\{\sigma_{\min}({\mathscr{M}}_{j}({\bm{\mathcal{M}}}))\right\\}$.
See remark after Theorem 2.

######  Theorem 2.

Let $\widehat{}\mathbf{s}^{(0)}$ be the initial clustering output by Algorithm
2. There exists some absolute constant $c,C_{1},C_{2},C_{3},C_{4}>0$ such that
if

$\displaystyle\Lambda_{\textsf{\tiny min}}\geq C_{1}(rK)^{1/2}d^{1/2}n^{1/4},$
(13)

and

$\displaystyle\Delta^{2}\geq{C_{2}\alpha^{-1}K^{2}}\left(\frac{dKr}{n}+1\right),$
(14)

we get, with probability at least $1-\exp(-c(n\wedge d))$, that

$n^{-1}\cdot h_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})\leq
C_{3}\frac{{K}}{\Delta^{2}}\left(\frac{dKr}{n}+1\right),$

and

$\displaystyle
n^{-1}\cdot\ell_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})\leq
C_{4}{\gamma^{2}K}\left(\frac{dKr}{n}+1\right),$

where $\gamma:=\max_{a\neq
b\in[K]}\left\|\mathbf{M}_{a}-\mathbf{M}_{b}\right\|_{\rm F}/\Delta$.

Theorem 2 suggests that Algorithm 2 delivers a consistent clustering if the
separation strength $\Delta^{2}\gg K(1+rdK/n)$. In terms of loss function
$\ell_{c}(\cdot)$, we have an additional dependence on $\gamma$, which relates
$\Delta$ to maximum separation strength. A similar condition is also casted in
the vector GMM (Lu and Zhou, 2016). As argued in Lu and Zhou (2016); Jin et
al. (2016), a distant cluster can cause local search to fail which indicates
the possibly unavoidable dependence on $\gamma$. Furthermore, Theorem 2
imposes a condition on the tensor signal strength $\Lambda_{\textsf{\tiny
min}}$, which is not needed in Theorem 1. Such an eigen-gap type condition is
prevalent in low-rank models (Zhang and Xia, 2018; Richard and Montanari,
2014; Levin et al., 2019; Xia, 2021; Lyu and Xia, 2022) as it determines
whether the population centers or their singular spaces are estimable by
polynomial-time algorithms, only in which case the low-rank structure can be
beneficial. Remarkably, $\Lambda_{\textsf{\tiny min}}$ also governs the
computational and statistical limit under LrMM as will be explained in Section
4.

Finally, by combining Theorem 2 and Theorem 1, the successes of Algorithm 1
and Algorithm 2 require the signal strength and separation strength conditions

$\Lambda_{\textsf{\tiny min}}\geq C_{1}(rK)^{1/2}d^{1/2}n^{1/4}$

and

$\quad\frac{\Delta^{2}}{\alpha^{-1}\gamma^{2}(\kappa_{0}^{2}\vee
Kr)Kr\left(\frac{dKr}{n}+1\right)}\rightarrow\infty$

To facilitate a clearer understanding of $\Lambda_{\textsf{\tiny min}}$, we
introduce the concept of individual signal strength denoted by $\lambda$. This
quantity, which is common in low-rank matrix literature, is defined as the
minimum value of the smallest singular value among $\mathbf{M}_{k}$’s, i.e.,

$\lambda:=\min_{k\in[K]}\sigma_{\textsf{\tiny min}}(\mathbf{M}_{k})$

Relation between tensor signal strength $\Lambda_{\textsf{\tiny min}}$ and
individual matrix signal strength $\lambda$. Define the condition number of
${\bm{\mathcal{M}}}$ in the mode-$j$ as
$\kappa_{j}:=\left\|{\mathscr{M}}_{j}({\bm{\mathcal{M}}})\right\|/\sigma_{\min}({\mathscr{M}}_{j}({\bm{\mathcal{M}}}))$
for $j=1,2$.

######  Lemma 1.

For $j\in\\{1,2\\}$,
$\sigma_{\min}({\mathscr{M}}_{j}({\bm{\mathcal{M}}}))\geq\kappa_{j}^{-1}(Kr)^{-1/2}\sqrt{n}\lambda$.

By Lemma 1, a sufficient condition for (13) to hold can be casted as
$\lambda\geq C_{0}(\kappa_{1}\vee\kappa_{2})rKd^{1/2}n^{-1/4}$. Recall that
$\kappa_{0}$ tells whether individual population center matrices are well-
conditioned. Here $\kappa_{1}$ ($\kappa_{2}$, resp.) measures the goodness of
alignment among the column (row, resp.) spaces of all population center
matrices. However, the exact relation between $\kappa_{1}$ and the column
spaces $\\{{\rm ColSpan}(\mathbf{U}_{k}^{\ast})\\}_{k=1}^{K}$ can be
intricate. The following lemma unfolds two special cases. Recall that
$r_{\mathbf{U}}$ and $r_{\mathbf{V}}$ are the ranks of
$\mathbf{U}=(\mathbf{U}_{1},\cdots,\mathbf{U}_{K})$ and
$\mathbf{V}=(\mathbf{V}_{1},\cdots,\mathbf{V}_{K})$, respectively, and
$\mathring{r}=\sum_{k=1}^{K}r_{k}$. Denote $\kappa(\mathbf{U})$ and
$\kappa(\mathbf{V})$ the condition numbers of $\mathbf{U}$ and $\mathbf{V}$,
respectively. The following indicates the connection between $\kappa_{j}$ and
$\kappa_{0}$.

######  Lemma 2.

Let ${\bm{\mathcal{M}}}$ admits low-rank decomposition (8). We have

$\displaystyle{\mathscr{M}}_{1}({\bm{\mathcal{M}}}){\mathscr{M}}_{1}^{\top}({\bm{\mathcal{M}}})=$
$\displaystyle\mathbf{U}\cdot{\rm
diag}\big{(}\\{n_{k}^{\ast}\mathbf{\Sigma}_{k}^{2}\\}_{k=1}^{K}\big{)}\cdot\mathbf{U}^{\top}$
$\displaystyle{\mathscr{M}}_{2}({\bm{\mathcal{M}}}){\mathscr{M}}_{2}^{\top}({\bm{\mathcal{M}}})=$
$\displaystyle\mathbf{V}\cdot{\rm
diag}\big{(}\\{n_{k}^{\ast}\mathbf{\Sigma}_{k}^{2}\\}_{k=1}^{K}\big{)}\cdot\mathbf{V}^{\top}$

and $\kappa_{1}\leq\kappa_{0}\kappa(\mathbf{U})\cdot(n^{\ast}_{\textsf{\tiny
max}}/n^{\ast}_{\textsf{\tiny min}})^{1/2}$ and
$\kappa_{2}\leq\kappa_{0}\kappa(\mathbf{V})\cdot(n^{\ast}_{\textsf{\tiny
max}}/n^{\ast}_{\textsf{\tiny min}})^{1/2}$ where $n^{\ast}_{\textsf{\tiny
min}}:=\min_{k}n_{k}^{\ast}$ and $n^{\ast}_{\textsf{\tiny
max}}:=\max_{k}n_{k}^{\ast}$. If $r_{\mathbf{U}}=r_{\mathbf{V}}=r_{1}$, i.e.,
all the population center matrices share the same singular space with
$\mathbf{M}_{1}$, we have
$\max\\{\kappa_{1},\kappa_{2}\\}\leq\kappa_{0}\cdot(K^{2}/\alpha)^{1/2};$ if
$r_{\mathbf{U}}=r_{\mathbf{V}}=\mathring{r}$ and $\mathbf{M}_{k}$ has mutually
orthogonal singular space, we have
$\max\\{\kappa_{1},\kappa_{2}\\}\leq\kappa_{0}\cdot(K/\alpha)^{1/2}$.

According to Lemma 2, the unfolded matrices
${\mathscr{M}}_{1}({\bm{\mathcal{M}}})$ and
${\mathscr{M}}_{2}({\bm{\mathcal{M}}})$ are well-conditioned if $\mathbf{U}$
and $\mathbf{V}$ are well-conditioned. Interestingly, this implies that our
tensor-based spectral initialization becomes more efficient when the
population center matrices $\mathbf{M}_{k}$’s have either perfectly aligned
singular spaces or nearly orthogonal singular spaces.

Discussions on tensor signal strength $\Lambda_{\textsf{\tiny min}}$.
Condition (13) reflects the computational difficulty under LrMM. This
intrinsic computational condition is likely attributed to the tensor method,
which is solely present in the initialization stage (Algorithm 2). Once well
initialized, the requirement for $\Lambda_{\min}$ vanishes in Theorem 1 for
lr-Lloyd (Algorithm 1). Such conditions are common in tensor problems Zhang
and Xia (2018); Auddy and Yuan (2022); Richard and Montanari (2014); Luo and
Zhang (2022). A more relevant work Lyu and Xia (2022) provides evidence
showing that no polynomial time can consistently estimate the population
centers even in the symmetric two-component LrMM if $\Lambda_{\textsf{\tiny
min}}=o(d^{1/2}n^{1/4})$. In Section 4, evidences are provided showing that
the same phenomenon exists for clustering, that is, if $\Lambda_{\textsf{\tiny
min}}=o(d^{1/2}n^{1/4})$, consistent clustering is impossible by any
polynomial time algorithms even when the separation strength $\Delta$ is much
stronger than the minimal condition (11).

Comparison with HOOI (Zhang and Xia, 2018) and the condition number of
${\mathscr{M}}_{3}({\bm{\mathcal{M}}})$.  Algorithm 2 looks similar to HOOI
(Zhang and Xia, 2018), which uses HOSVD for mode-wise spectral initialization
and applies power iterations to further improve the estimates of singular
spaces. Indeed, (13) is analogous to the signal strength condition for HOOI
therein to succeed. However, the mode-wise HOSVD and subsequent power
iterations both require a lower bound on $\sigma_{\textsf{\tiny
min}}\big{(}{\mathscr{M}}_{k}({\bm{\mathcal{M}}})\big{)},k=1,2,3$. While our
Theorem 2 also requires a lower bound on $\sigma_{\textsf{\tiny
min}}\big{(}{\mathscr{M}}_{1}({\bm{\mathcal{M}}})\big{)}$ and
$\sigma_{\textsf{\tiny
min}}\big{(}{\mathscr{M}}_{2}({\bm{\mathcal{M}}})\big{)}$, we emphasize that a
similar lower bound on $\sigma_{\textsf{\tiny
min}}\big{(}{\mathscr{M}}_{3}({\bm{\mathcal{M}}})\big{)}$ is too strong and
trivialize the whole problem. To see this, just notice via definition that
$\Delta\geq\sigma_{\textsf{\tiny
min}}\big{(}{\mathscr{M}}_{3}({\bm{\mathcal{M}}})\big{)}/2$.

Comparison with Han et al. (2022a). A tensor block model was proposed by Han
et al. (2022a), which can be regarded as an extension of the stochastic block
model. They developed the high-order Lloyd’s algorithm (HLloyd) with spectral
initialization. The two works differ drastically from several aspects. From
the algorithmic perspective, HLloyd doesn’t require low-rank approximation at
all since it explores block structure rather than low-rank structure. The
membership matrix in Han et al. (2022a) (analogous to $\mathbf{U}_{k}$ in this
paper) lies in the space $\\{0,1\\}^{d_{k}\times r_{k}}$, which is more
informative owing to its discrete structure. Clearly, block model is just a
special case of low-rank model and HLloyd is inapplicable to our LrMM. On the
technical front, HLolyd updates the block means simply by the sample average
which admits an explicit and clean representation form. In sharp contrast, the
analysis for lr-Llyod is much more challenging due to the implicit and
complicated form of the updated cluster centers
$\widehat{\mathbf{M}}_{k}^{(t)}$ defined in (6), which calls for more advanced
tools.

###  3.3 Minimax lower bound

Theorem 1 has shown that the low-rank Lloyd’s algorithm achieves the
asymptotical clustering error rate $\exp(-\Delta^{2}/8)$. In this section, a
matching minimax lower bound is derived showing that the aforesaid rate is
indeed optimal in the minimax sense. A lower bound under GMM has been
established by Lu and Zhou (2016). We follow the arguments in Gao et al.
(2018) to establish the minimax lower bound for LrMM. Observe that the error
rate only depends on the separation strength $\Delta$ implying that the
dimension $d_{1},d_{2}$ and ranks $r_{k}$’s play a less important role here.

Define the following parameter space for the population center matrices and
arrangements of latent labels:

$\displaystyle\Omega_{\Delta}\equiv\Omega(\Delta,d_{1},d_{2},n,K,\alpha):=\Big{\\{}(\\{$
$\displaystyle\mathbf{M}_{k}\\}_{k=1}^{K},\mathbf{s}):~{}\mathbf{M}_{k}\in\mathbb{R}^{d_{1}\times
d_{2}},\text{rank}(\mathbf{M}_{k})=r_{k},\mathbf{s}\in[K]^{n},$
$\displaystyle\min_{k\in[K]}|\\{i\in[n]:s_{i}=k\\}|\geq\alpha n/K,\min_{a\neq
b}\left\|\mathbf{M}_{a}-\mathbf{M}_{b}\right\|_{\rm F}\geq\Delta\Big{\\}}$

For notation simplicity, we omit its dependence on the ranks $r_{k}$’s.

######  Theorem 3.

Let $\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ satisfy LrMM (3) with
$(\\{\mathbf{M}_{k}\\}_{k=1}^{K},\mathbf{s}^{\ast})\in\Omega_{\Delta}$.
Suppose $\\{\mathbf{E}_{i}\\}_{i=1}^{n}$ has i.i.d ${\cal N}(0,\sigma^{2})$
entries. If
${\Delta^{2}}/\left(\sigma^{2}\log(K/\alpha)\right)\rightarrow\infty$ as
$n\to\infty$, we have

$\inf_{\widehat{}\mathbf{s}}\sup_{(\\{\mathbf{M}_{k}\\}_{k=1}^{n},\mathbf{s}^{\ast})\in\Omega_{\Delta}}\
\mathbb{E}\frac{h_{\textsf{c}}(\widehat{}\mathbf{s},\mathbf{s}^{\ast})}{n}\geq\exp\left(-(1+o(1))\frac{\Delta^{2}}{8\sigma^{2}}\right)$

where $\underset{\widehat{}\mathbf{s}}{\inf}$ is taken over all clustering
algorithms.

Compared to Theorem 1 and Theorem 2, the minimax lower bound is established
only requiring a separation strength $\Delta^{2}\gg 1$ assuming
$K/\alpha=O(1)$. Theorem 3 holds for any signal strength and the infimum is
taking over all possible clustering algorithms without considering their
computational feasibility. Here, an algorithm is said computationally feasible
if it is computable within a polynomial time complexity in terms of $n$ and
$d_{1},d_{2}$.

##  4 Computational Barriers

We now turn to the computational hardness of LrMM. For simplicity, we set
$\alpha,K,r\asymp 1$ throughout this section. Our signal strength condition
(13) in initialization requires a lower bound $\Lambda_{\textsf{\tiny
min}}\gtrsim d^{1/2}n^{1/4}$. The purpose of this section is to provide
evidences on its necessity to guarantee computationally feasible clustering
algorithms. Our evidence is built on the low-degree likelihood ratio framework
for hypothesis testing proposed by Kunisky et al. (2019); Hopkins (2018),
which has delivered convincing evidences justifying the computational hardness
under sparse GMM (Löffler et al., 2020) and for sparse PCA (Ding et al.,
2019).

Suppose that, given i.i.d. observations
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$, one is interested in the computational
and statistical limit in distinguishing two hypothesis $\mathbb{Q}_{n}$ and
$\mathbb{P}_{n}$, i.e,

$H_{0}^{(n)}:\mathbf{X}_{1}\sim\mathbb{Q}_{n}\quad{\rm versus}\quad
H_{1}^{(n)}:\mathbf{X}_{1}\sim\mathbb{P}_{n}$ (15)

The above two hypotheses are said statistically indistinguishable if no test
can have both type I and type II error probabilities vanishing asymptotically.
The famous Neyman-Pearson lemma tells us that the likelihood ratio test based
on
$L_{n}({\bm{\mathcal{X}}}):=d\mathbb{P}_{n}/d\mathbb{Q}_{n}(\mathbf{X}_{1},\cdots,\mathbf{X}_{n})$
has a preferable power and is uniformly most powerful under some scenarios. A
well recognized fact is that $\mathbb{Q}_{n}$ and $\mathbb{P}_{n}$ are
statistically indistinguishable if the quantity
$\|L_{n}\|^{2}:=\mathbb{E}_{\mathbb{Q}_{n}}[L_{n}({\bm{\mathcal{X}}})^{2}]$
remains bounded as $n\to\infty$. See Kunisky et al. (2019) for a simple proof.

While the asymptotic magnitude of $\|L_{n}\|^{2}$ is informative for
understanding the statistical limit of testing (15), it does not directly
reflect the computational limit of testing (15). Towards that end, the low-
degree likelihood ratio framework seeks a polynomial approximation of
$L_{n}({\bm{\mathcal{X}}})$ and investigates the magnitude of the resultant
approximation. More exactly, let $L_{n}^{\leq D}({\bm{\mathcal{X}}})$ be the
orthogonal projection of $L_{n}({\bm{\mathcal{X}}})$ onto the linear space
spanned by polynomials $\mathbb{R}^{d_{1}\times d_{2}\times
n}\mapsto\mathbb{R}$ of degrees at most $D$. Similarly, define $\|L_{n}^{\leq
D}\|^{2}:=\mathbb{E}_{\mathbb{Q}_{n}}[L_{n}^{\leq
D}({\bm{\mathcal{X}}})^{2}]$. Kunisky et al. (2019) conjectures that the
asymptotic magnitude of $\|L_{n}^{\leq D}\|^{2}$ reflects the computational
hardness of testing the hypothesis (15). More formally, their conjecture,
slightly adapted for our purpose, can be written as follows. It has been
introduced in Lyu and Xia (2022). Here, a test $\phi_{n}(\cdot)$ taking value
$1$ means rejecting the null hypothesis and takes value $0$ if the null
hypothesis is not rejected. Thus
$\mathbb{E}_{\mathbb{Q}_{n}}[\phi_{n}({\bm{\mathcal{X}}})]$ and
$\mathbb{E}_{\mathbb{P}_{n}}[1-\phi_{n}({\bm{\mathcal{X}}})]$ stands for
type-I and type-II error, respectively.

######  Conjecture 1 (Lyu and Xia (2022)).

If there exists $\epsilon>0$ and $D=D_{n}\geq(\log nd)^{1+\epsilon}$ for which
$\left\|L_{n}^{\leq D}\right\|=1+o(1)$ as $n\to\infty$, then there is no
polynomial-time test $\phi_{n}:\mathbb{R}^{d_{1}\times d_{2}\times
n}\mapsto\\{0,1\\}$ such that the sum of type-I error and type-II error
probabilities

$\mathbb{E}_{\mathbb{Q}_{n}}[\phi_{n}({\bm{\mathcal{X}}})]+\mathbb{E}_{\mathbb{P}_{n}}[1-\phi_{n}({\bm{\mathcal{X}}})]\rightarrow
0\quad\text{as}\quad n\rightarrow\infty$

Based on this conjecture, Kunisky et al. (2019) reproduces the sharp phase
transitions for the spiked Wigner matrix model and the widely-believed
statistical-to-computational gap in tensor PCA, and Lyu and Xia (2022)
develops a computational hardness theory for estimating the population low-
rank matrices under LrMM.

Note that a specific hypothesis $\mathbb{P}_{n}$ is necessary to apply
Conjecture 1 and investigate the computational barriers in clustering for
LrMM. Towards that end, we consider a symmetric two-component LrMM as in Lyu
and Xia (2022). It is a special case of model (3) with $K=2$, $r_{1}=r_{2}=1$,
$\mathbf{M}_{1}=n^{-1/2}\Lambda_{\textsf{\tiny
min}}\mathbf{u}\mathbf{v}^{\top}$ and
$\mathbf{M}_{2}=-\mathbf{M}_{1}=-n^{-1/2}\Lambda_{\textsf{\tiny
min}}\mathbf{u}\mathbf{v}^{\top}$. Here $\mathbf{u}\in\mathbb{R}^{d_{1}}$ and
$\mathbf{v}\in\mathbb{R}^{d_{2}}$ have unit norms. In this case, the tensor
signal strength is $\Lambda_{\textsf{\tiny min}}>0$. Moreover, the individual
signal strength is $\lambda=n^{-1/2}\Lambda_{\textsf{\tiny min}}$ and
separation strength is $\Delta=2n^{-1/2}\Lambda_{\textsf{\tiny min}}$, i.e.,
the two quantities are at the same order. Then the observations can be re-
written as

$\mathbf{X}_{i}=s_{i}^{\ast}(n^{-1/2}\Lambda_{\textsf{\tiny
min}}\mathbf{u}\mathbf{v}^{\top})+\mathbf{E}_{i},\quad\forall i=1,\cdots,n,$
(16)

where $s_{i}^{\ast}=1$ if $\mathbf{X}_{i}$ is sampled from ${\cal
N}(\mathbf{M}_{1},\mathbf{I}_{d_{1}}\otimes\mathbf{I}_{d_{2}})$ and
$s_{i}^{\ast}=-1$ if $\mathbf{X}_{i}$ is sampled from ${\cal
N}(\mathbf{M}_{2},\mathbf{I}_{d_{1}}\otimes\mathbf{I}_{d_{2}})$. Note that the
rank-one model (16) is no more difficult than the general K-component case but
it suffices for our purpose. The null hypothesis $\mathbb{Q}_{n}$ corresponds
to the case $\Lambda_{\textsf{\tiny min}}=0$, i.e., all observations are pure
noise. Clearly, the difficulty level of distinguishing $\mathbb{Q}_{n}$ and
$\mathbb{P}_{n}$ is characterized by signal strength $\Lambda_{\textsf{\tiny
min}}$ in eq. (16). Conjecture 1 requires the calculation of $\|L_{n}^{\leq
D}\|^{2}$, which is extremely difficult for generally fixed singular vectors
$\mathbf{u},\mathbf{v}$ and deterministic latent labels $\mathbf{s}^{\ast}$. A
prior distribution simplifies the calculation. Finally, our null and
alternative hypothesis are formally defined as follows.

######  Definition 1 (Null and alternative hypothesis).

  * •

Under $\mathbb{Q}_{n}$, we observe $n$ matrices
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ generated i.i.d. from (16) with
$\Lambda_{\textsf{\tiny min}}=0$. Equivalently, it means that each
$\mathbf{X}_{i}$ has i.i.d. standard normal entries.

  * •

Under $\mathbb{P}_{n}:=\mathbb{P}_{n}^{\Lambda_{\textsf{\tiny min}}}$, we
observe $n$ matrices $\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ generated i.i.d.
from (16) with $\Lambda_{\textsf{\tiny min}}>0$, and moreover, each coordinate
of $\mathbf{u}$ and $\mathbf{v}$ independently uniformly take values from
$\\{\pm d_{1}^{-1/2}\\}$ and $\\{\pm d_{2}^{-1/2}\\}$, respectively, and the
entries of $\mathbf{s}^{\ast}$ are independent Rademacher random variables,
i.e., taking $\pm 1$ with equal probabilities.

######  Theorem 4.

Consider $\mathbb{Q}_{n}$ and $\mathbb{P}_{n}$ in Definition 1. If
$\Lambda_{\textsf{\tiny min}}=o\left(d^{1/2}n^{1/4}\right)$ as $n\to\infty$,
then $\left\|L_{n}^{\leq D}\right\|=1+o(1)$.

The proof of Theorem 4 can be found in Lyu and Xia (2022). If Conjecture 1 is
true, Theorem 4 implies that $\mathbb{Q}_{n}$ and
$\mathbb{P}_{n}^{\Lambda_{\textsf{\tiny min}}}$ are statistically
indistinguishable by polynomial-time algorithms as long as the signal strength
$\Lambda_{\textsf{\tiny min}}=o\left(d^{1/2}n^{1/4}\right)$. We now establish
the connection of testing the hypothesis to the clustering problem under two-
component symmetric LrMM (16).

For any fixed $\Lambda_{\textsf{\tiny min}}>0$, define the parameter space of
interest by

$\displaystyle\widetilde{\Omega}_{\Lambda_{\textsf{\tiny min}}}$
$\displaystyle\equiv\widetilde{\Omega}(\Lambda_{\textsf{\tiny
min}},d_{1},d_{2},n)$
$\displaystyle=\Big{\\{}(\mathbf{M},\mathbf{s}):~{}\mathbf{M}=n^{-1/2}\Lambda^{\prime}_{\textsf{\tiny
min}}\mathbf{u}\mathbf{v}^{\top},\mathbf{u}\in\mathbb{R}^{d_{1}},\mathbf{v}\in\mathbb{R}^{d_{2}},\mathbf{s}\in\\{\pm
1\\}^{n},|{\bf 1}^{\top}\mathbf{s}|\leq n/2,\Lambda^{\prime}_{\textsf{\tiny
min}}\geq\Lambda_{\textsf{\tiny min}}\Big{\\}}$

By Chernoff bound, with probability at least $1-e^{-c_{0}n}$ where $c_{0}>0$
is an absolute constant, the i.i.d. observations
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}$ generated by
$\mathbb{P}_{n}^{\Lambda_{\textsf{\tiny min}}}$ satisfy the rank-one LrMM (16)
with parameters
$(\mathbf{M},\mathbf{s})\in\tilde{\Omega}_{\Lambda_{\textsf{\tiny min}}}$. The
following theorem tells that if consistent clustering is possible for LrMM, so
is for distinguishing the hypothesis in Definition 1.

######  Theorem 5.

Suppose there exists a clustering algorithm
$\widehat{}\mathbf{s}_{\textsf{comp}}:\mathbb{R}^{d_{1}\times d_{2}\times
n}\mapsto\\{\pm 1\\}^{n}$ for LrMM (16) with runtime $poly(n,d)$ that is
consistent under the sequence of signal strength
$\big{\\{}\Lambda_{\textsf{\tiny min}}^{(n)}\big{\\}}_{n\geq 1}$ in the sense
that there exists a sequence $\\{(\delta_{n},\zeta_{n})\\}_{n\geq 1}\to 0$
such that for all large $n$,

$\sup_{(\mathbf{M},\mathbf{s}^{\ast})\in\widetilde{\Omega}_{\Lambda_{\textsf{\tiny
min}}^{(n)}}}\mathbb{P}\left(n^{-1}\cdot
h_{\textsf{c}}(\widehat{}\mathbf{s}_{\textsf{comp}},\mathbf{s}^{\ast})>\delta_{n}\right)\leq\zeta_{n}$
(17)

If the signal strength satisfies $\Lambda_{\textsf{\tiny min}}^{(n)}\geq
C_{0}(1+\epsilon^{-2})^{1/2}d^{1/2}$ with some absolute constant $C_{0}>0$ and
$\epsilon\in(0,1)$, then there exists a test $\phi_{n}:\mathbb{R}^{d_{1}\times
d_{2}\times n}\mapsto\\{0,1\\}$ with runtime $poly(n,d)$ that consistently
distinguishes $\mathbb{P}_{n}^{\Lambda_{\textsf{\tiny min}}^{(n)}}$ from
$\mathbb{Q}_{n}$ so that

$\mathbb{E}_{Q_{n}}[\phi_{n}({\bm{\mathcal{X}}})]+\sup_{((1-\epsilon)\mathbf{M},\mathbf{s}^{\ast})\in\widetilde{\Omega}_{\Lambda_{\textsf{\tiny
min}}^{(n)}}}\mathbb{E}_{(\mathbf{M},\mathbf{s}^{\ast})}[1-\phi_{n}({\bm{\mathcal{X}}})]\to
0,\quad{\rm as}\ n,d\to\infty.$

Essentially, Theorem 5 only needs a signal strength $\Lambda_{\textsf{\tiny
min}}\gg d^{1/2}$ to successfully reduce a polynomial-time clustering
algorithm to a polynomial-time hypothesis test. Based on Conjecture 1, a
combination of Theorem 4 and Theorem 5 implies the following result, whose
proof is straightfoward and hence omitted.

######  Corollary 1.

Suppose Conjecture 1 holds for $\mathbb{Q}_{n}$ and $\mathbb{P}_{n}$ in
Definition 1. If the signal strength $\Lambda_{\textsf{\tiny
min}}^{(n)}=o(d^{1/2}n^{1/4})$, then for any polynomial-time clustering
algorithm $\widehat{}\mathbf{s}_{\textsf{comp}}$, there exist absolute
constants $\delta,\zeta>0$ such that

$\sup_{(\mathbf{M},\mathbf{s}^{\ast})\in\widetilde{\Omega}_{\Lambda_{\textsf{\tiny
min}}^{(n)}}}\mathbb{P}\left(n^{-1}\cdot
h_{\textsf{c}}(\widehat{}\mathbf{s}_{\textsf{comp}},\mathbf{s}^{\ast})>\delta\right)\geq\zeta$

as $n\to\infty$.

It is worth pointing out that even though the signal strength
$\Lambda_{\textsf{\tiny min}}=o(d^{1/2}n^{1/4})$, the separation strength
$\Delta=2n^{-1/2}\Lambda_{\textsf{\tiny min}}$ can still be much larger than
$d^{1/2}n^{-1/2}$ that is required by Theorem 1. This suggests that if signal
strength is not strong, consistent clustering by polynomial-time algorithms is
still impossible even though the separation strength is very strong.

##  5 Relaxing the Signal Strength Condition

Our main theorem in Section 3 imposes a strong signal strength condition on
all the population center matrices, i.e., $\Lambda_{\textsf{\tiny min}}$ is
lower bounded by $\Omega(d^{1/2}n^{1/4})$, or equivalently, $\lambda$ is lower
bounded by $\Omega(d^{1/2}n^{-1/4})$. While evidences in Section 4 show that
this condition might be necessary for the two-component symmetric case if only
polynomial-time algorithms are sought, this condition appears flawed in the
general asymmetric case. This section aims to relax the signal strength
condition in the sense that one population center matrix is allowed to be
arbitrarily smaller (in spectral norm) than $d^{1/2}n^{-1/4}$, in which case
(13) might fail.

To simplify the narrative, we focus on the two-component LrMM, i.e., $K=2$ in
model (3), whose population center matrices are denoted by $\mathbf{M}_{1}$
and $\mathbf{M}_{2}$, respectively. However, it is straightforward to extend
our discussion to the general case. For $K=2$, it is more intuitive and
convenient to express everything in terms of individual signal strength
$\mathbf{M}_{1}$ and $\mathbf{M}_{2}$ instead of the tensor signal strength
$\Lambda_{\min}$, even though they are equivalent444Alternatively, we can
impose condition on $\min_{j=1,2}\\{\sigma_{\textsf{\tiny
min}}({\bm{\mathcal{M}}}_{j})\\}$, where $[{\bm{\mathcal{M}}}_{j}]_{\cdot\cdot
i}=\mathbb{I}(s_{i}^{\ast}=1)\mathbf{M}_{j}$.. Without loss of generality, we
assume that $\|\mathbf{M}_{1}\|_{\rm F}$ is large so that reliable estimation
is possible, and that $\|\mathbf{M}_{2}\|_{\rm F}$ is small so that reliable
estimation is impossible. The following assumption is made to clarify this
further.

######  Assumption 2.

There exists a small constant $c>0$ such that

$\sigma_{1}(\mathbf{M}_{2})\leq
c\alpha^{-1/2}\left(\sqrt{\frac{d}{n}}+\kappa_{0}^{-1}\right),$

and

$\frac{\sigma^{2}_{r_{1}}(\mathbf{M}_{1})}{\alpha^{-1}(\kappa_{0}^{2}\vee
r_{1})\left(\frac{d}{n}+1\right)}\rightarrow\infty$

where $\kappa_{0}$, with slight abuse of notation, is the condition number of
$\mathbf{M}_{1}$.

If $\kappa_{0},\alpha=O(1)$, Assumption 2 can be recasted as
$\sigma_{r_{1}}(\mathbf{M}_{1})\gg d^{1/2}n^{-1/2}+1$ and
$\sigma_{1}(\mathbf{M}_{2})\leq c(d^{1/2}n^{-1/2}+1)$. Note that Assumption 2
puts no lower bound on $\sigma_{1}(\mathbf{M}_{2})$. In the extreme case,
$\sigma_{1}(\mathbf{M}_{2})$ is allowed to be zero and consistent estimation
of $\mathbf{M}_{2}$ is unavailable even if the true labels are revealed.
Assumption 2 already implies that $\Delta\gg\left(d^{1/2}n^{-1/2}+1\right)$ if
the ranks $r_{1},r_{2}$ are both upper bounded by $O(1)$, matching the
separation condition (11) in Theorem 1. Intuitively, although clustering shall
becomes easier as the constant $c$ in Assumption 2 decreases, this cannot be
verified by Theorem 1 where the signal strength condition (13) fails.

Under Assumption 2, it is generally pointless to compute the center matrix
$\widehat{}\mathbf{M}_{2}$ by SVD in Lloyd’s algorithm since $\mathbf{M}_{2}$
cannot be reliably estimated. Moreover, the SVD procedure complicates the
subsequent theoretical analysis of Lloyd’s algorithm. Instead of estimating
$\mathbf{M}_{2}$ via SVD, we opt to a trivial estimate by setting
$\widehat{}\mathbf{M}_{2}^{(t)}={\bf 0}$. The detailed steps are enumerated in
Algorithm 3, whose theoretical performance is guaranteed by Theorem 6.555We
remark that the low-rankness assumption for $\mathbf{M}_{2}$ in Theorem 6 is
not essential, which can be dropped by instead requiring
${\sqrt{r_{1}}\sigma_{r_{1}}(\mathbf{M}_{1})}/\left\|\mathbf{M}_{2}\right\|_{\rm
F}\to\infty$.

Algorithm 3 Low-rank Lloyd’s Algorithm under Relaxed SNR Assumption 2 (rlr-
Lloyd)

Input: Observations:
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\in\mathbb{R}^{d_{1}\times d_{2}}$ where
$\mathbf{X}_{i}=\mathbf{M}_{s_{i}^{\ast}}+\mathbf{E}_{i}$ and
$s_{i}^{\ast}\in\\{1,2\\}$, initial estimate $\widehat{}\mathbf{s}^{(0)}$,
ranks $r_{1},r_{2}$.

for $t=1,\ldots,T$ do

For each $k=1,2$:  $\widehat{}\mathbf{M}_{k}^{(t)}\leftarrow\textrm{best
rank-}r_{k}\textrm{ approximation of
}\bar{\mathbf{X}}_{k}(\widehat{}\mathbf{s}^{(t-1)}):=\frac{\sum_{i=1}^{n}\mathbb{I}\left(\widehat{s}_{i}^{(t-1)}=k\right)\mathbf{X}_{i}}{\sum_{i=1}^{n}\mathbb{I}\left(\widehat{s}_{i}^{(t-1)}=k\right)}$

Set $\widehat{}\mathbf{M}_{2}^{(t)}\leftarrow{\bf 0}$ if
$\sigma_{1}(\widehat{}\mathbf{M}_{2}^{(t)})<\sigma_{1}(\widehat{}\mathbf{M}_{1}^{(t)})$;
or set
$\widehat{}\mathbf{M}_{1}^{(t)}\leftarrow\widehat{}\mathbf{M}_{2}^{(t)}$,
$\widehat{}\mathbf{M}_{2}^{(t)}\leftarrow{\bf 0}$ if
$\sigma_{1}(\widehat{}\mathbf{M}_{2}^{(t)})>\sigma_{1}(\widehat{}\mathbf{M}_{1}^{(t)})$.

Re-label by setting, for each $i\in[n]$:
$\widehat{s}_{i}^{(t)}\leftarrow\underset{k\in[2]}{\arg\min}\
\|\mathbf{X}_{i}-\widehat{}\mathbf{M}_{k}^{(t)}\|_{\rm F}^{2}$

end for

Output: $\widehat{}\mathbf{s}=\widehat{}\mathbf{s}^{(T)}$

######  Theorem 6.

Suppose Assumption 2 holds and $d\geq C_{0}\log K$ for some absolute constant
$C_{0}>0$. Assume $\widehat{}\mathbf{s}^{(0)}$ satisfies

$n^{-1}\cdot\ell_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})=o\left(\frac{\alpha}{\kappa_{0}^{2}}\Delta^{2}\right)$
(18)

Furthermore, if
$\sqrt{\frac{r_{1}}{r_{2}}}\cdot\frac{\sigma_{r_{1}}(\mathbf{M}_{1})}{\sigma_{1}(\mathbf{M}_{2})}\to\infty$,
then we have

$n^{-1}\cdot
h_{\textsf{c}}(\widehat{}\mathbf{s}^{(t-1)},\mathbf{s}^{\ast})\leq\exp\left(-\big{(}1-o(1)\big{)}\frac{\Delta^{2}}{8}\right)+\frac{1}{2^{t}}$

with probability at least $1-\exp(-\Delta)-\exp\big{(}-c_{0}d\big{)}$ for a
small absolute constant $c_{0}>0$.

To ensure a consistent $\widehat{}\mathbf{s}^{(0)}$ satisfying (18), we use a
modified version of the tensor initialization discussed in Section 3.2. The
original spectral initialization can be mis-leading if a rank $r_{\mathbf{U}}$
larger than $r_{1}$ is adopted. For our purpose, only the top-$r_{1}$ singular
vectors are taken during spectral initialization, i.e., effort is made only
for estimating $\mathbf{M}_{1}$ whose left and right singular vectors are
denoted by $\mathbf{U}_{1}$ and $\mathbf{V}_{1}$, respectively. See Algorithm
4 for further algorithmic details and Theorem 7 for theoretical guarantees.

Algorithm 4 Tensor-based Spectral Initialization Under Relaxed SNR Assumption
(rTS-Init)

Input: observations:
$\mathbf{X}_{1},\cdots,\mathbf{X}_{n}\in\mathbb{R}^{d_{1}\times d_{2}}$ where
$\mathbf{X}_{i}=\mathbf{M}_{s_{i}^{\ast}}+\mathbf{E}_{i}$ and
$s_{i}^{\ast}\in\\{1,2\\}$; or a tensor
${\bm{\mathcal{X}}}\in\mathbb{R}^{d_{1}\times d_{2}\times n}$ by concatenating
the matrix observations slice by slice, ranks $r_{1}$.

Spectral initialization:

  * 1.

Obtain the estimated singular vectors $\widehat{}\mathbf{U}_{1}$ and
$\widehat{}\mathbf{V}_{1}$ by applying HOSVD to the tensor
${\bm{\mathcal{X}}}$ in mode-1 and mode-2 matricizations with rank $r_{1}$.

  * 2.

Project ${\bm{\mathcal{X}}}$ onto the column space of
$\widehat{}\mathbf{U}_{1}$ and $\widehat{}\mathbf{V}_{1}$ by
$\widehat{\bm{\mathcal{G}}}:={\bm{\mathcal{X}}}\times_{1}\widehat{}\mathbf{U}_{1}\widehat{}\mathbf{U}_{1}^{\top}\times_{2}\widehat{}\mathbf{V}_{1}\widehat{}\mathbf{V}_{1}^{\top}$

  * 3.

Apply K-means on the rows of
$\widehat{}\mathbf{G}:={\mathscr{M}}_{3}(\widehat{\bm{\mathcal{G}}})\in\mathbb{R}^{n\times
d_{1}d_{2}}$ and obtain the initial clustering by

$(\widehat{}\mathbf{s}^{(0)},\\{\widehat{}\mathbf{M}_{1}^{(0)},\widehat{}\mathbf{M}_{2}^{(0)}\\}):=\underset{\mathbf{s}\in[2]^{n};\mathbf{M}_{1},\mathbf{M}_{2}\in\mathbb{R}^{d_{1}\times
d_{2}}}{\arg\min}\
\sum_{i=1}^{n}\big{\|}[\widehat{}\mathbf{G}]_{i\cdot}-vec(\mathbf{M}_{s_{i}})\big{\|}^{2}$

Output: $\widehat{}\mathbf{s}^{(0)}$

######  Theorem 7.

Let $\widehat{}\mathbf{s}^{(0)}$ be the initial clustering output by Algorithm
4. Suppose there exists constant $c,C_{0}>0$ and large constant $C>1$ such
that $n/\kappa_{0}^{4}\geq C$,

$\sigma_{r_{1}}(\mathbf{M}_{1})\geq
C\alpha^{-1/2}\frac{d^{1/2}}{n^{1/4}},\quad\sigma_{1}(\mathbf{M}_{2})\leq
C^{-1}\kappa_{0}^{-1}\frac{d^{1/2}}{n^{1/4}},$

then we get, with probability at least $1-\exp(-cd)$, that

$n^{-1}\cdot
h_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})\leq\frac{C_{0}}{\Delta^{2}}\left(\frac{dr_{1}}{n}+1\right)\quad\text{and}\quad
n^{-1}\cdot\ell_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})\leq{C_{0}}\left(\frac{dr_{1}}{n}+1\right).$

Furthermore, if $n/\kappa_{0}^{4}\rightarrow\infty$ and
$\alpha\Delta^{2}/\kappa_{0}^{2}\rightarrow\infty$, with probability at least
$1-\exp(-cd)$ we have that

$n^{-1}\cdot
h_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})=o\left(\frac{\alpha}{\kappa_{0}^{2}}\right)\quad\text{and}\quad
n^{-1}\cdot\ell_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})=o\left(\frac{\alpha}{\kappa_{0}^{2}}\Delta^{2}\right).$

Theorem 7 serves as a counterpart of Theorem 2, with the distinction that we
express the conditions in terms of $\sigma_{r_{1}}(\mathbf{M}_{1})$ and
$\sigma_{1}(\mathbf{M}_{2})$. Notably, the threshold $d^{1/2}n^{-1/4}$
illuminates the disparity between statistical and computational aspects in the
presence of low-rankness structure as discussed in Section 4. We emphasize
that the gap arises solely due to the initialization procedure similar to the
case in Section 3.2. Within our framework, Assumption 2 together with good
initializer $\widehat{}\mathbf{s}^{(0)}$ suffices to guarantee the statistical
optimality of Algorithm 3 under relaxed a signal strength condition and
minimal requirement on the separation strength $\Delta$.

##  6 Clustering versus Estimation

Lyu and Xia (2022) investigated the minimax optimal estimation of latent low-
rank matrices under two-component symmetric LrMM, which revealed multiple
phase transitions and a statistical-to-computational gap. In this section,
together with Theorem 1 and 2, we discuss the differences between estimation
and clustering.

###  6.1 Example where clustering is more challenging

For simplicity, we consider the rank-one symmetric two-component LrMM (16)
with $d_{1}=d_{2}=d$, where the separation strength $\Delta\asymp
n^{-1/2}\Lambda_{\min}$ and individual signal strength
$\lambda=n^{-1/2}\Lambda_{\min}$ coincides up to a constant factor. To make
comparison, in this section we consider $\Lambda_{\textsf{\tiny min}}$ instead
of $\lambda$. The minimax rate of estimating $\mathbf{M}$ (up to a sign flip),
established in Lyu and Xia (2022), is

$\inf_{\widehat{}\mathbf{M}}\sup_{(\mathbf{M},\mathbf{s}^{\ast})\in\widetilde{\Omega}_{\Lambda_{\textsf{\tiny
min}}}}\mathbb{E}\min_{\eta=\pm
1}\left\|\widehat{}\mathbf{M}-\eta\mathbf{M}\right\|_{\rm
F}\asymp\min\left\\{d^{1/2}{\Lambda^{-1}_{\textsf{\tiny
min}}}+d^{1/2}n^{-1/2},n^{-1/2}\Lambda_{\textsf{\tiny min}}\right\\}$ (19)

The above rate is achievable by the computationally NP-hard maximum likelihood
estimator with almost no constraint on signal strength and by a
computationally fast spectral-aggregation estimator under the regime of strong
signal strength $\Lambda_{\textsf{\tiny min}}\gtrsim d^{1/2}n^{1/4}$. For a
fair comparison, we focus on this computationally feasible regime. The phase
transitions under this regime can be summarized as in Table 2.

Sample size | Signal strength | Minimax optimal estimation error
---|---|---
$d^{2}\lesssim n$ | $d^{1/2}n^{1/4}\lesssim\Lambda_{\textsf{\tiny min}}\lesssim n^{1/2}$ |  $\frac{\sqrt{d}}{\Lambda_{\textsf{\tiny min}}}$
$\Lambda_{\textsf{\tiny min}}\gtrsim n^{1/2}$ |  $\sqrt{\frac{d}{n}}$
$d^{2}\gg n$ | $\Lambda_{\textsf{\tiny min}}\gtrsim d^{1/2}n^{1/4}$ |  $\sqrt{\frac{d}{n}}$

Table 2: Phase transition in minimax optimal estimation for two-component
symmetric LrMM under the regime of strong signal strength
$\Lambda_{\textsf{\tiny min}}\gtrsim d^{1/2}n^{1/4}$. See (19) and Lyu and Xia
(2022) for more details.

Without loss of generality, we assume the dimension $d\to\infty$ as
$n\to\infty$. The case $d^{2}\gg n$ is referred to as the high-dimensional
setting, and $d^{2}\lesssim n$ is called the low-dimensional setting. An
estimator $\widehat{}\mathbf{M}$ is said strongly consistent if the relative
estimation error $\|\widehat{}\mathbf{M}-\mathbf{M}\|_{\rm
F}\|\mathbf{M}\|_{\rm F}^{-1}$ approaches to zero in expectation as
$n\to\infty$. Table 2 tells that strongly consistent estimation $\mathbf{M}$
is always achievable as long as the signal strength is greater
$d^{1/2}n^{1/4}$. A particularly interesting regime is
$d^{1/2}n^{1/4}\lesssim\Lambda_{\textsf{\tiny min}}\lesssim n^{1/2}$. For
instance, when $d^{2}=o(n)$, $\mathbf{M}$ can still be consistently estimated
even when the signal strength $\Lambda_{\textsf{\tiny min}}\to 0$ as
$n\to\infty$.

It is certainly not the case for clustering. Besides consistent clustering
(see definition in Theorem 5), we say a clustering algorithm is weakly
efficient if it can beat a random guess, but the mis-clustering error rate
does not vanish as $n\to\infty$. When $d^{2}=o(n)$, Theorem 3 dictates that
even weakly efficient clustering is impossible, i.e.,
$\exp(-\Lambda_{\textsf{\tiny min}}^{2}/(2n))$ is at least $1/2$, if
$\Lambda_{\textsf{\tiny min}}\leq c_{0}n^{1/2}$ for some absolute constant
$c_{0}>0$. However, the spectral aggregation estimator (Lyu and Xia, 2022) can
still consistently estimate the population center matrix $\mathbf{M}$ in the
aforesaid scenario. Moreover, by Theorem 1, consistent clustering even
requires $\Lambda_{\textsf{\tiny min}}/n^{1/2}\to\infty$, which is much more
stringent than that required by (strongly) consistent estimation.

The differences of phase transitions in estimation and clustering are
enumerated in Table 3. Basically, strongly consistent estimation is always
possible as long as $\Lambda_{\textsf{\tiny min}}\gtrsim d^{1/2}n^{1/4}$. In
contrast, weakly efficient clustering is possible only when
$\Lambda_{\textsf{\tiny min}}\gtrsim n^{1/2}+d^{1/2}n^{1/4}$, and consistent
clustering is possible only when $\Lambda_{\textsf{\tiny min}}\gtrsim
d^{1/2}n^{1/4}$ and meanwhile $\Lambda_{\textsf{\tiny min}}\gg n^{1/2}$. Note
that the gap between estimation and clustering is present only under the low-
dimensional setting $n\gtrsim d^{2}$. The gap vanishes under the high-
dimensional setting $d^{2}\gg n$, in which case the signal strength condition
$\Lambda_{\textsf{\tiny min}}\gtrsim d^{1/2}n^{1/4}$ already implies
$\Lambda_{\textsf{\tiny min}}\gg n^{1/2}$.

Sample size | Signal strength | Consistent estimation | Weakly efficient clustering | Consistent clustering
---|---|---|---|---
$d^{2}\lesssim n$ | $d^{1/2}{n^{1/4}}\lesssim\Lambda_{\textsf{\tiny min}}\lesssim n^{1/2}$ | Possible | Impossible | Impossible
$n^{1/2}\lesssim\Lambda_{\textsf{\tiny min}}\lesssim n^{1/2}$ | Possible | Possible | Impossible
$\Lambda_{\textsf{\tiny min}}\gg n^{1/2}$ | Possible | Possible | Possible
$d^{2}\gg n$ | $\Lambda_{\textsf{\tiny min}}\gtrsim d^{1/2}{n^{1/4}}$ | Possible | Possible | Possible

Table 3: The differences of phase transitions in estimation and clustering for
two-component symmetric LrMM under the regime of strong signal strength
$\Lambda_{\textsf{\tiny min}}\gtrsim d^{1/2}n^{1/4}$. Here $d^{2}\gg n$ is
referred to as the high-dimensional setting, and $d^{2}\lesssim n$ as the low-
dimensional setting.

We collect these facts to convince that, at least for the two-component
symmetric LrMM (16), clustering is intrinsically more challenging than
estimation. The same phenomenon also arises in GMM. See, e.g., Wu and Zhou
(2019).

###  6.2 Example where estimation is more challenging

While, generally, clustering is recognized as being more challenging than
estimation, there are examples where clustering is easier than estimation.
Similarly as in Section 5, consider the two-component LrMM with population
center matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$ so that

$\sigma_{r_{1}}(\mathbf{M}_{1})\geq
C_{1}\Big{(}1+\frac{d^{1/2}}{n^{1/2}}+\frac{d^{1/2}}{n^{1/4}}\Big{)}\quad{\rm
and}\quad\sigma_{1}(\mathbf{M}_{2})\leq
C_{1}^{-1}\cdot\frac{d^{1/2}}{n^{1/2}}$

where $C_{1}>0$ is a large constant and, for simplicity, we assume
$\kappa_{0},\alpha,r_{1},r_{2}=O(1)$. Observe that

$\sqrt{\frac{r_{1}}{r_{2}}}\cdot\frac{\sigma_{r_{1}}(\mathbf{M}_{1})}{\sigma_{1}(\mathbf{M}_{2})}\geq\begin{cases}C_{1}^{2}n^{1/4},&\textrm{
if }n\leq d^{2};\\\ C_{1}^{2}(n/d)^{1/2},&\textrm{ if
}n>d^{2};\end{cases}\quad\to\infty,\quad{\rm as}\ n\to\infty$

Moreover,

$\Delta:=\|\mathbf{M}_{1}-\mathbf{M}_{2}\|_{\rm F}\gtrsim
C_{1}\Big{(}1+\frac{d^{1/2}}{n^{1/4}}\Big{)}\to\infty$

if the constant $C_{1}>0$ diverges to infinity. Therefore, by Theorem 6, if
$C_{1}\to\infty$, our Algorithm 3 consistently cluster all observations.

However, consistent estimation of the population center matrices is more
challenging. Even if all the latent labels are correctly identified,
estimation of $\mathbf{M}_{2}$ is still impossible because of its weak signal
strength. Indeed, the low-rank approximation to

$\bar{\mathbf{X}}_{2}(\mathbf{s}^{\ast}):=\frac{1}{n_{2}^{\ast}}\sum_{i=1}^{n}\mathbb{I}\left(s_{i}^{\ast}=2\right)\mathbf{X}_{i}$

achieves the error rate (in expectation) $O(d^{1/2}n^{-1/2})$ and the relative
error rate (in expectation) diverges to infinity as $C_{1}\to\infty$.
Similarly, the trivial estimate by a zero matrix attains the relative error
rate $1$ that never vanishes as $n\to\infty$. Consequently, a strongly
consistent estimate of $\mathbf{M}_{2}$ becomes impossible.

##  7 Discussions

###  7.1 Estimation of $r_{\mathbf{U}}$, $r_{\mathbf{V}}$, $K$ and $r_{k}$’s

Our tensor-based spectral initialization method requires an input of ranks
$r_{\mathbf{U}}$, $r_{\mathbf{V}}$ and the number of clusters $K$, which are
usually unknown in practice. Under the decomposition (9), they constitute the
Tucker ranks of tensor ${\bm{\mathcal{M}}}$. Several approaches are available
to estimate the Tucker ranks for tensor PCA model. One typical approach (Jing
et al., 2021; Cai et al., 2022) is to check the scree plots (Cattell, 1966) of
${\mathscr{M}}_{1}({\bm{\mathcal{X}}})$,
${\mathscr{M}}_{2}({\bm{\mathcal{X}}})$ and
${\mathscr{M}}_{3}({\bm{\mathcal{X}}})$, respectively. Under a suitable signal
strength condition as in Theorem 2, the scree plots of
${\mathscr{M}}_{1}({\bm{\mathcal{X}}})$ and
${\mathscr{M}}_{2}({\bm{\mathcal{X}}})$ shall serve a reliable estimate of
$r_{\mathbf{U}}$ and $r_{\mathbf{V}}$, respectively. However, we note that it
is statistically more efficient to estimate $K$ by, instead, taking the scree
plot of
${\mathscr{M}}_{3}({\bm{\mathcal{X}}}\times_{1}\widehat{}\mathbf{U}^{\top}\times_{2}\widehat{}\mathbf{V}^{\top})$,
where $\widehat{}\mathbf{U}$ and $\widehat{}\mathbf{V}$ are obtained in step 1
of Algorithm 2. This additional spectral projection promotes further noise
reduction as in Algorithm 2. After obtaining $r_{\mathbf{U}}$,
$r_{\mathbf{V}}$ and $K$, an initial clustering $\widehat{}\mathbf{s}^{(0)}$
can be attained by apply Algorithm 2. Similarly, we then estimate the rank
$r_{k}$ by the scree plot of the sample average of matrix observations whose
initial labels are $k$. It provides a valid estimate as long as the initial
clustering is sufficiently good. The aforementioned approach works nicely in
real-world data applications. See Section 8 for more details.

###  7.2 Matrix observation with categorical entries

Oftentimes, the matrix observations consist of categorical entries. For
instance, the Malaria parasite gene networks (see Section 8.2.3) have binary
entries (Bernoulli distribution); the 4D-scanning transmission electron
microscopy (Han et al., 2022b) produces count-type entries (Poisson
distribution). Our algorithms are still applicable and deliver appealing
performance on, e.g., Malaria parasite gene networks dataset. Unfortunately,
our theory can not directly cover those cases, although the noise are still
sub-Gaussian. Without loss of generality, let us consider multi-layer binary
networks and assume $\mathbf{X}_{i}$ has Bernoulli entries. Then the entries
of $\mathbf{X}_{i}$ have an equal variance only when they have the same
expectation, reducing the network to a trivial Erdős-Rényi graph.
Nevertheless, equal noise variance is crucial to establish Theorem 2.
Moreover, the techniques for proving Theorem 1 are likely sub-optimal since
the sub-Gaussian constant $\sigma_{\textsf{sg}}$ is usually not sharp enough
to characterize a Bernoulli random variable. We leave this to future works.

##  8 Numerical Experiments and Real Data Applications

###  8.1 Numerical Experiments

This section presents the empirical performance of lr-Lloyd’s algorithm
(Algorithm 1) and its relaxed variant under weak SNR (Algorithm 3) referred to
as the rlr-Lloyd’s algorithm. Specifically, we focus on the algorithmic
convergence and final clustering error.

In the first simulation setting S1, we fix the dimension $d_{1}=d_{2}=50$ and
sample size $n=200$. The latent labels $s_{i}^{\ast}$ are generated i.i.d.
from the model (2) with equal mixing probabilities, i.e., $\pi_{k}=1/K$. All
the presented results in S1 are based on the average of $30$ independent
trials. We test the convergence of Algorithm 1 under both Gaussian (S1-1) and
Bernoulli (S1-2) noise.

In S1-1, we set $K=2$, $r_{1}=r_{2}=2$ and standard Gaussian noise. The
population center matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$ are generated
in the following manner. For each $k=1,2$, we independently generate a
$d_{1}\times d_{2}$ matrix with i.i.d. standard Gaussian entries and extract
its top-$2$ left and right singular vectors as $\mathbf{U}_{k}$ and
$\mathbf{V}_{k}$, respectively. The singular values are manually set as
$\mathbf{\Sigma}_{k}=\text{diag}\\{1.2\lambda,\lambda\\}$ for some fix
$\lambda>0$. Then the population center matrices are constructed as
$\mathbf{M}_{k}=\mathbf{U}_{k}\mathbf{\Sigma}_{k}\mathbf{V}_{k}^{\top}$. Our
experiment tries four levels of signal strength
$\lambda\in\\{1.9,2.1,2.3,2.5\\}$. For each $\lambda$, the population center
matrices are generated as above and the separation strength is recorded. The
corresponding separation strength are $\Delta\in\\{4.22,4.66,5.11,5.45\\}$. At
each level of signal strength, the observations
$\\{\mathbf{X}_{i}:i=1,\cdots,200\\}$ are independently drawn from (4) with
the obtained center matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$. Here we
focus on the convergence behavior of Lloyd’s iterations of Algorithm 1, and
thus a warm initial clustering $\widehat{}\mathbf{s}^{(0)}$ is provided before
hand. The same initial clustering is used for all simulations and the initial
clustering error is
$n^{-1}h_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})=0.45$,
i.e., slightly better than a random guess. Convergence of Algorithm 1 under
four levels of signal strength (or, correspondingly, separation strength) is
displayed in the left plot of Figure 1. The decreasing of log of clustering
error is linear in first few iterations, as expected by our Theorem 1. The
algorithm converges fast and the final clustering error is reflected by the
separation strength $\Delta$. It is worth pointing out that Figure 1(a) also
shows that Algorithm 1 converges faster when $\Delta$ becomes larger. While
this cannot be directly concluded from Theorem 1, it can be easily verified by
checking the proof.

In S1-2, we test the effectiveness of Algorithm 1 under non-Gaussian and
non-i.i.d. noise. In particular, we consider the mixture multi-layer
stochastic block model (MMSBM) introduced in Jing et al. (2021) 666We
emphasize that our Theorem 1 is not directly applicable to MMSBM due to
non-i.i.d. noise.. We set the number of clusters $K=3$. For each $k=1,2,3$,
the $k$-th SBM is associated with a connection probability matrix
$\mathbf{B}_{k}\in[0,1]^{K\times K}$ and a membership matrix
$\mathbf{Z}_{k}\in\\{0,1\\}^{d\times K}$, which are set as
$\mathbf{B}_{k}:=\bar{p}_{k}\cdot\mathbf{I}_{K}+\bar{p}_{k}/2\cdot(\mathbf{1}_{K}\mathbf{1}_{K}^{\top}-\mathbf{I}_{K})$
with $\bar{p}_{k}=\bar{p}\cdot k/K$ and
$\mathbf{Z}_{k}(i,:)=\mathbf{e}_{s^{*}_{i}}$, respectively. Thus each SBM has
three cluster of nodes and the population center matrices are
$\mathbf{M}_{k}=\mathbf{Z}_{k}\mathbf{B}_{k}\mathbf{Z}_{k}^{\top}\in[0,1]^{d\times
d}$. Conditioned on the latent label $\mathbf{s}_{i}^{\ast}$, the $i$-th
observation $\mathbf{X}_{i}$ is sampled from
SBM($\mathbf{Z}_{s_{i}^{\ast}},\mathbf{B}_{s_{i}^{\ast}}$), namely,
$\mathbf{X}_{i}(j_{1},j_{2})\sim\text{Bernoulli}(\mathbf{M}_{s^{*}_{i}}(j_{1},j_{2}))$
and $\mathbf{X}_{i}(j_{2},j_{1})=\mathbf{X}_{i}(j_{1},j_{2})$ for $1\leq
j_{1}<j_{2}\leq d$. Note that $\mathbf{X}_{i}$ is symmetric because the
network is undirected. We manually set the diagonal entries of
$\mathbf{X}_{i}$ to zeros so that no self-loop is allowed in the observed
network. Clearly, the entry-wise variances of $\mathbf{X}_{i}$ are not
necessarily equal. Under the above MMSBM, the signal strength and separation
strength are characterized by sparsity level $\bar{p}$. Four sparsity levels
$\bar{p}\in\\{0.05,0.08,0.10,0.15\\}$ are studied so that the corresponding
separation strength are $\Delta\in\\{0.75,1.19,1.46,2.15\\}$. Similarly, a
fixed good initial clustering $\widehat{}\mathbf{s}^{(0)}$ is used for all
simulations and the initial clustering error is
$n^{-1}h_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})=0.3$.
Convergence behavior of Algorithm 1 is displayed in the right plot of Figure
1. Still, Lloyd’s iterations converges fast and the final clustering error is
decided by the separation strength $\Delta$.

(a) Simulation S1-1: Log of clustering error ($K=2$) with $\Delta$ varying
under Gaussian noise.

(b) Simulation S1-2 Log of clustering error ($K=3$) with $\Delta$ varying
under Bernoulli noise (MMSBM).

Figure 1: (Convergence behavior of Algorithm 1) Log of clustering error with
$\Delta$ varying under two scenarios: LrMM with Gaussian noise and MMSBM with
Bernoulli noise.

In the second simulation setting S2, we aim to compare the final clustering
error of vanilla Lloyd’s algorithm and our low-rank Lloyd’s algorithm. The
dimensions are varied at two cases $d_{1}=d_{2}\in\\{50,100\\}$, sample size
is set as $n\in\\{100,200\\}$, number of clusters $K=2$ and ranks
$r_{1}=r_{2}=3$. The latent labels are generated as in S1. For each $d_{1}$
and $n$, the simulation is repeated for 100 times and their average clustering
error rate is reported.

In S2-1, the population center matrices $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$
are constructed such that they share identical singular spaces. More exactly,
we extract singular vectors $\mathbf{U}_{1}$, $\mathbf{V}_{1}$ and singular
value matrix $\mathbf{\Sigma}_{1}$ as is done in S1-1. Then the population
center matrices are set as
$\mathbf{M}_{1}=\mathbf{U}_{1}\mathbf{\Sigma}_{1}\mathbf{V}_{1}^{\top}$ and
$\mathbf{M}_{2}=\mathbf{U}_{1}(\mathbf{\Sigma}_{1}+\text{diag}\\{\Delta/3,\Delta/3,\Delta/3\\})\mathbf{V}_{1}^{\top}$.
Here the signal strength is fixed at $\lambda=10$ and the separation parameter
is chosen from $\Delta\in\\{1,5,10\\}$. The final clustering error and its
standard error by four methods are reported in the upper half of Table 4.
Noted that the initialization of “vec-Lloyd” in Lu and Zhou (2016) is attained
by spectral clustering on ${\mathscr{M}}_{3}({\bm{\mathcal{X}}})$. We observe
that the clustering errors of four methods all decrease as $\Delta$ increases.
However, lr-Lloyd initialized by Algorithm 2 achieves a much smaller
clustering error compared with other methods. This is due to the fact that our
proposed tensor-based spectral initialization is capable to capture the low-
rank signal whereas both spectral clustering and naive K-means on
${\mathscr{M}}_{3}({\bm{\mathcal{X}}})$ ignores the low-rank structure in the
other two modes of ${\bm{\mathcal{M}}}$. As a result, all the other three
methods perform almost the same under current setting. Lastly, the bold-font
column in Table 4 confirms Theorem 1 in that the clustering error achieved by
TS-init initialized lr-Lloyd algorithm is only determined by $\Delta$
regardless of the dimension $d_{1},d_{2}$ or the sample size $n$.

In S2-2, the singular vectors of $\mathbf{M}_{1}$ and $\mathbf{M}_{2}$ are
generated exactly the same as in S1-1. The singular values of $\mathbf{M}_{1}$
and $\mathbf{M}_{2}$ are set as $\mathbf{\Sigma}_{1}={\rm
diag}(1.2\lambda,1.1\lambda,\lambda)$ and $\mathbf{\Sigma}_{2}={\rm
diag}(0.36,0.33,0.30)$, respectively. Then $\sigma_{\textsf{\tiny
min}}(\mathbf{M}_{1})=\lambda$ and $\sigma_{1}(\mathbf{M}_{2})=0.36$. Here
$\lambda$ is varied at $\\{1.9,2.2,2.5\\}$ for the case $d_{1}=d_{2}=50$ and
$\\{2.7,3.0,3.3\\}$ for the case $d_{1}=d_{2}=100$. Consequently, the signal
strength of $\mathbf{M}_{2}$ is much smaller than $\mathbf{M}_{1}$ that
corresponds to the weak SNR setting in Section 5, and we test the performance
of the relaxed lr-Lloyd’s algorithm (Algorithm 3). The results are reported in
the lower half of Table 4. Clearly, rlr-Lloyd’s algorithm outperforms the
vanilla Lloyd’s algorithm (i.e., the vectorized version). In certain cases,
the vanilla Lloyd’s algorithm merely beats a random guess whereas the rlr-
Lloyd’s algorithm almost achieves zero clustering error. We also observe that
rlr-Lloyd’s algorithm still performs nicely if initialized by K-means on
${\mathscr{M}}_{3}({\bm{\mathcal{X}}})$.

Setting | $d_{1}=d_{2}$ | $n$ | $\lambda$ | $\Delta$ |  | vec-Lloyd
---
(Lu and Zhou, 2016)
| lr-Lloyd initialized by
---
TS-Init (Algorithm 2)
| vec-Lloyd initialized by
---
K-means on ${\mathscr{M}}_{3}({\bm{\mathcal{X}}})$
| lr-Lloyd initialized by
---
K-means on ${\mathscr{M}}_{3}({\bm{\mathcal{X}}})$
S2-1 | 50 | 100 | 10 | 1 | 0.461 (0.032) | 0.401 (0.058) | 0.462 (0.030) | 0.459 (0.031)
10 | 5 | 0.459 (0.033) | 0.163 (0.039) | 0.456 (0.033) | 0.452 (0.034)
10 | 10 | 0.458 (0.034) | 0.066 (0.025) | 0.441 (0.047) | 0.433 (0.054)
200 | 10 | 1 | 0.475 (0.019) | 0.398 (0.056) | 0.469 (0.025) | 0.466 (0.025)
10 | 5 | 0.473 (0.021) | 0.152 (0.027) | 0.462 (0.027) | 0.450 (0.039)
10 | 10 | 0.471 (0.022) | 0.063 (0.016) | 0.437 (0.041) | 0.380 (0.082)
100 | 100 | 10 | 1 | 0.461 (0.028) | 0.391 (0.069) | 0.460 (0.033) | 0.461 (0.033)
10 | 5 | 0.461 (0.029) | 0.157 (0.054) | 0.455 (0.036) | 0.455 (0.036)
10 | 10 | 0.460 (0.029) | 0.063 (0.026) | 0.458 (0.034) | 0.456 (0.034)
200 | 10 | 1 | 0.468 (0.023) | 0.390 (0.064) | 0.469 (0.023) | 0.467 (0.023)
10 | 5 | 0.468 (0.024) | 0.147 (0.028) | 0.469 (0.022) | 0.465 (0.026)
10 | 10 | 0.467 (0.024) | 0.062 (0.017) | 0.459 (0.030) | 0.451 (0.037)
Setting | $d_{1}=d_{2}$ | $n$ | $\sigma_{\textsf{\tiny min}}(\mathbf{M}_{1})$ | $\Delta$ |  | vec-Lloyd
---
(Lu and Zhou, 2016)
| rlr-Lloyd
---
(Algorithm 3)
| vec-Lloyd initialized by
---
K-means on ${\mathscr{M}}_{3}({\bm{\mathcal{X}}})$
| rlr-Lloyd initialized by
---
K-means on ${\mathscr{M}}_{3}({\bm{\mathcal{X}}})$
S2-2 | 50 | 100 | 1.9 | 3.68 | 0.434 (0.052) | 0.314 (0.138) | 0.418 (0.066) | 0.327 (0.129)
2.2 | 4.24 | 0.424 (0.061) | 0.134 (0.125) | 0.385 (0.079) | 0.152 (0.138)
2.5 | 4.81 | 0.417 (0.068) | 0.041 (0.051) | 0.309 (0.103) | 0.055 (0.091)
200 | 1.9 | 3.68 | 0.433 (0.052) | 0.070 (0.020) | 0.380 (0.070) | 0.072 (0.046)
2.2 | 4.24 | 0.431 (0.054) | 0.057 (0.018) | 0.351 (0.077) | 0.059 (0.048)
2.5 | 4.81 | 0.424 (0.057) | 0.035 (0.015) | 0.268 (0.088) | 0.033 (0.014)
100 | 100 | 2.7 | 5.19 | 0.422 (0.056) | 0.300 (0.169) | 0.416 (0.057) | 0.301 (0.164)
3 | 5.76 | 0.421 (0.059) | 0.131 (0.164) | 0.390 (0.077) | 0.176 (0.181)
3.3 | 6.33 | 0.426 (0.053) | 0.067 (0.139) | 0.347 (0.086) | 0.065 (0.130)
200 | 2.7 | 5.19 | 0.442 (0.040) | 0.019 (0.010) | 0.395 (0.071) | 0.022 (0.037)
3 | 5.76 | 0.443 (0.041) | 0.008 (0.006) | 0.301 (0.089) | 0.008 (0.007)
3.3 | 6.33 | 0.440 (0.043) | 0.003 (0.004) | 0.190 (0.069) | 0.003 (0.004)

Table 4: Clustering error of lr-Lloyd (Algorithm 1) and rlr-Lloyd (Algorithm
3) compared with vanilla Lloyd’s algorithm (Lu and Zhou, 2016) on vectorized
data (vec-Lloyd). The number in brackets represents the standard error over
$100$ trials.

###  8.2 Real Data Applications

We now demonstrate the merits of our proposed low-rank Lloyd’s (lr-Lloyd)
algorithm on several real-world datasets and compare with existing methods.

####  8.2.1 BHL dataset

The BHL (brain, heart and lung ) dataset777The dataset is publicly available
at https://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS1083., which had been
analyzed in Mai et al. (2021), consists of $d_{1}=1124$ gene expression
profiles of $n=27$ brain, heart, or lung tissues. Each tissue is measured
repeatedly for $d_{2}=4$ times and hence the $i$th sample can be constructed
as $\mathbf{X}_{i}\in\mathbb{R}^{1124\times 4}$ for $i=1,\cdots,27$. Our aim
is to correctly identify those $\mathbf{X}_{i}$’s belonging to the same type
of tissue, i.e., $K=3$. We apply Algorithm 1 together with an initial
clustering $\widehat{}\mathbf{s}^{(0)}$ obtained by Algorithm 2 with
$r_{\mathbf{U}}=r_{\mathbf{V}}=1$. These ranks are chosen based on the scree
plots of ${\mathscr{M}}_{1}({\bm{\mathcal{X}}})$ and
${\mathscr{M}}_{2}({\bm{\mathcal{X}}})$. The final clustering error attained
by lr-Lloyd’s algorithm is $n^{-1}\cdot
h_{\textsf{c}}(\widehat{}\mathbf{s},\mathbf{s}^{\ast})=0.03704$. As shown in
Table 5, our lr-Lloyd’s algorithm performs the best among all the
competitors888Note that all results except lr-Lloyd are directly borrowed from
Mai et al. (2021), which use $\mathbf{X}_{i}$’s after dimension reduction to a
size of either $20\times 4$ or $30\times 4$, and we only report the better one
here. that are reported in Mai et al. (2021).

| lr-Lloyd | DEEM | K-means | SKM | DTC | TBM | EM | AFPF
---|---|---|---|---|---|---|---|---
Clustering error | 3.70 | 7.41 | 11.11 | 11.11 | 18.52 | 11.11 | 11.11 | 11.11
Table 5: Clustering error on BHL dataset. SKM: sparse K-means (Witten and
Tibshirani, 2010); DTC: dynamic tensor clustering (Sun and Li, 2019); TBM:
tensor block model (TBM) (Wang and Zeng, 2019); EM: standard EM implemented in
Mai et al. (2021); AFPF: adaptive pairwise fusion penalized clustering (Guo et
al., 2010).

The improvement can be attributed to two reasons. First, DEEM in Mai et al.
(2021) is designed based on EM algorithm targeted at Gaussian probability
distribution, and hence they need to first perform multiple Kolmogorov-Smirnov
tests to drop the columns not following Gaussian distribution, which might
lead to potential information loss. In sharp contrast, their procedure is not
necessary for our method, as the low-rank Lloyd’s algorithm allows for sub-
Gaussian noise. Secondly, our algorithm is more suitable for the specific
structure of the data. Particularly, the population center matrices are
expected to be rank-one as the columns of $\mathbf{X}_{i}$ represent repeated
measurements for the same sample. However, such planted structure is under-
exploited in Mai et al. (2021) and others.

####  8.2.2 EEG dataset

The EEG dataset999The dataset is publicly available at
https://archive.ics.uci.edu/ml/datasets/EEG+Database. has been extensively
studied by various statistical models (Li et al., 2010; Zhou and Li, 2014; Hu
et al., 2020; Huang et al., 2022). The goal is to inspect EEG correlations of
genetic predisposition to alcoholism. The data contains measurements which
were sampled at $d_{1}=256$ Hz for $1$ second, from $d_{2}=64$ electrodes
placed on each scalp of $n=122$ subjects. Each subject, either being alcoholic
or not, completed 120 trials under different stimuli. More detailed
description of the dataset can be found in Zhang et al. (1995). For our
application, we average all the trials for each subject under single stimulus
condition (S1) and two matched stimuli condition (S2), respectively, and
construct the data tensor as
${\bm{\mathcal{X}}}^{(S_{1})}\in\mathbb{R}^{256\times 64\times 122}$ (or
${\bm{\mathcal{X}}}^{(S_{2})}\in\mathbb{R}^{256\times 64\times 122}$) after
standardization. Thus each subject is associated with a $256\times 64$ matrix,
and we aim to cluster these subjects into $K=2$ groups, corresponding to
alcholic group and control group. We apply rlr-Lloyd’s algorithm (Algorithm 3)
with $r_{\mathbf{U}}=r_{\mathbf{V}}=3$ and $r_{1}=2,r_{2}=1$. Here
$r_{\mathbf{U}}$ and $r_{\mathbf{V}}$ are selected by the scree plot of
${\mathscr{M}}_{1}({\bm{\mathcal{X}}})$ and
${\mathscr{M}}_{2}({\bm{\mathcal{X}}})$, and $r_{1}$ and $r_{2}$ are tuned by
interpreting the final outcomes. The clustering error of our method and
competitors are shown in Table 6. It is worth pointing out that our task of
clustering is generally more challenging than classification, which has been
investigated on the EEG dataset (Li et al., 2010; Zhou and Li, 2014; Hu et
al., 2020; Huang et al., 2022). Those classification approaches often achieve
lower classification error rates. As a faithful comparison, our rlr-Lloyd’s
algorithm enjoys a superior performance to its competitors in terms of
clustering error rate and time complexity.

Surprisingly, we note that the original lr-Lloyd’s algorithm (Algorithm 1 \+
Algorithm 2) would not deliver a satisfactory result on this dataset. It can
be partially explained by Figure 2, which displays the average of all trials
under S2 for two groups. It is readily seen that the average matrix of control
group is comparatively close to pure noise, and hence the relaxed version of
lr-Lloyd’s algorithm can work reasonably well in this scenario.

| rlr-Lloyd | vec-Lloyd | SKM | DTC | TBM
---|---|---|---|---|---
S1 | 39.34 | 42.62 | 44.26 | 45.08 | 43.44
S2 | 28.69 | 35.25 | 36.07 | 39.34 | 35.25
Table 6: Clustering error of EEG dataset under S1 and S2. Note that the
methods vec-Lloyd and SKM (Witten and Tibshirani, 2010) refer to directly
applying Lloyd’s algorithm and sparse K-means on vectorized data, i.e., on
rows of ${\mathscr{M}}_{3}({\bm{\mathcal{X}}}^{(S_{1})})$ or
${\mathscr{M}}_{3}({\bm{\mathcal{X}}}^{(S_{2})})$, whereas DTC(Sun and Li,
2019) and TBM (Wang and Zeng, 2019) are both tensor-based clustering methods.

Figure 2: EEG dataset: average of matrix observations for alcoholic group
(left) and control group (right) under S2.

####  8.2.3 Malaria parasite genes networks dataset

We then consider the var genes networks of the human malaria parasite
Plasmodium falciparum constructed by Larremore et al. (2013) via mapping $n=9$
highly variable regions (HVRs) to a multi-layer network.

(a) HVR1.

(b) HVR2.

(c) HVR3.

(d) HVR4.

(e) HVR5.

(f) HVR6.

(g) HVR7.

(h) HVR8.

(i) HVR9.

Figure 3: Malaria parasite genes networks dataset: 9 highly variable regions
(HVRs) represented by their adjacency matrices (Jing et al., 2021)

Following the practice in Jing et al. (2021), we focus on $d_{1}=d_{2}=212$
common nodes appearing on all $9$ layers and obtain a multi-layer network
adjacency tensor ${\bm{\mathcal{X}}}\in\\{0,1\\}^{212\times 212\times 9}$ with
each layer being the associated adjacency matrix. Unfortunately, the method in
Larremore et al. (2013) needs to discard $3$ out of $9$ HVRs due to their
extreme sparse structures, referring to region $\\{2,3,4\\}$ in Figure 3. This
later had been remedied by the tensor-decomposition-based method TWIST in Jing
et al. (2021). In term of clustering all layers, we expect our algorithm would
have a comparable performance in contrast with the results in Jing et al.
(2021). Specifically, Jing et al. (2021) obtain a hierarchical structure with
$6$ clusters of all layers by repeatedly clustering the embedding vectors.
Following their practice, by setting
$(r_{\mathbf{U}},r_{\mathbf{V}},K)=(15,15,6)$, we apply Algorithm 2 on
${\bm{\mathcal{X}}}$, and find that the $9$ HVRs fall in to the following
clusters: $\\{1\\},\\{2,3,4,5\\},\\{6\\},\\{7\\},\\{8\\},\\{9\\}$. The result
is exactly the same as that in Jing et al. (2021) but our method avoid
repeated clustering. We remark that our tensor-based spectral initialization
already produces a good initial clustering on this dataset, and thus further
low-rank Lloyd’s iterations seem unnecessary. In sharp contrast, it would lead
to unsatisfactory result if we directly apply K-means with $K=6$ on the
embedding matrix obtained by TWIST. This further demonstrates the validity and
flexibility of our proposed lr-Lloyd’s algorithm.

####  8.2.4 UN comtrade trade flow networks dataset

In the last example, we consider the international commodity trade flow data
in $2019$ in terms of countries/regions and different types of commodities,
collected by Lyu et al. (2021) from UN comtrade Database101010The dataset is
publicly available at https://comtrade.un.org.. Following the data processing
procedure in Lyu et al. (2021), we pick out top $d_{1}=d_{2}=48$
countries/regions ranked by exports and obtain a weighted adjacency tensor
$\widetilde{\bm{\mathcal{X}}}\in\mathbb{R}^{48\times 48\times 97}$, where
$n=97$ layers represent different categories of commodities111111The
categories are based on 2-digit HS code in https://www.foreign-
trade.com/reference/hscode.htm.. The entry
$\widetilde{\bm{\mathcal{X}}}(i_{1},i_{2},i_{3})$ indicates the amount of
exports from country $i_{1}$ to country $i_{2}$ in terms of commodity type
$i_{3}$. To have a comparable magnitude across different entries, our data
tensor is obtained after transformation
${\bm{\mathcal{X}}}=\log(\widetilde{\bm{\mathcal{X}}}+1)$. We emphasize that
in Lyu et al. (2021) the edges of ${\bm{\mathcal{X}}}$ have to be further
converted to binary under their framework, which might cause undesirable
information loss. We apply Algorithm 1 that is initialized by Algorithm 2 with
parameters $(r_{\mathbf{U}},r_{\mathbf{V}},K)=(3,3,2)$ and
$(r_{1},r_{2})=(2,2)$. These choices produce most interpretable result as
summarized in Table 7. It is intriguing to notice that cluster 1 mainly
consists of products of low durability including animal & vegetable products
and part of foodstuffs, whereas cluster 2 contains most industrial products
that might indicate a trend of global trading. These findings are consistent
with Lyu et al. (2021).

Commodity cluster 1 | Commodity cluster 2
---|---
01-05 Animal & Animal Products (100%) | 15 Vegetable Products (13.73%)
06-14 Vegetable Products (86.27%) | 19-22 Foodstuffs (60.82%)
16-18, 23-24 Foodstuffs (39.18%) | 25,27 Mineral Products (86.68%)
26 Mineral Products (13.32%) | 28-30,32-35,38 Chemicals & Allied Industries (96.46%)
31,36-37 Chemicals & Allied Industries (3.54%) | 39-40 Plastics / Rubbers (100%)
41,43 Raw Hides, Skins, Leather, & Furs (23.01%) | 42 Raw Hides, Skins, Leather, & Furs (76.99%)
45-47 Wood & Wood Products (15.13%) | 44,48-49 Wood & Wood Products (84.87%)
50-55,57-58,60 Textiles (23.40%) | 56,59,61-63 Textiles (65.97%)
65-67 Footwear / Headgear (17.45%) | 64 Footwear / Headgear (82.55%)
75,78-81 Metals (6.44%) | 68-71 Stone / Glass (100%)
86,89 Transportation (5.50%) | 72-74,76,82-83 Metals (93.56%)
91-93,97 Miscellaneous (8.19%) | 84-85 Machinery / Electrical (100%)
| 87-88 Transportation (94.50%)
| 90,94-96,99 Miscellaneous (91.81%)

Table 7: Clustering result of UN comtrade network. The number in brackets is
the percentage of the amount of exports in the corresponding type of
commodity.

## References

  * Abbe et al. (2020) Emmanuel Abbe, Jianqing Fan, and Kaizheng Wang.  An $\ell_{p}$ theory of pca and spectral clustering.  _arXiv preprint arXiv:2006.14062_ , 2020.
  * Arroyo et al. (2021) Jesús Arroyo, Avanti Athreya, Joshua Cape, Guodong Chen, Carey E Priebe, and Joshua T Vogelstein.  Inference for multiple heterogeneous networks with a common invariant subspace.  _Journal of Machine Learning Research_ , 22(142):1–49, 2021.
  * Athreya et al. (2017) Avanti Athreya, Donniell E Fishkind, Minh Tang, Carey E Priebe, Youngser Park, Joshua T Vogelstein, Keith Levin, Vince Lyzinski, and Yichen Qin.  Statistical inference on random dot product graphs: a survey.  _The Journal of Machine Learning Research_ , 18(1):8393–8484, 2017.
  * Auddy and Yuan (2022) Arnab Auddy and Ming Yuan.  On estimating rank-one spiked tensors in the presence of heavy tailed errors.  _IEEE Transactions on Information Theory_ , 68(12):8053–8075, 2022.
  * Balakrishnan et al. (2017) Sivaraman Balakrishnan, Martin J Wainwright, and Bin Yu.  Statistical guarantees for the em algorithm: From population to sample-based analysis.  _The Annals of Statistics_ , 45(1):77–120, 2017\.
  * Cai et al. (2021) Biao Cai, Jingfei Zhang, and Will Wei Sun.  Jointly modeling and clustering tensors in high dimensions.  _arXiv preprint arXiv:2104.07773_ , 2021.
  * Cai et al. (2022) Jian-Feng Cai, Jingyang Li, and Dong Xia.  Generalized low-rank plus sparse tensor estimation by fast riemannian optimization.  _Journal of the American Statistical Association_ , (just-accepted):1–39, 2022.
  * Cai and Zhang (2018) T Tony Cai and Anru Zhang.  Rate-optimal perturbation bounds for singular subspaces with applications to high-dimensional statistics.  _The Annals of Statistics_ , 46(1):60–89, 2018\.
  * Cattell (1966) Raymond B Cattell.  The scree test for the number of factors.  _Multivariate behavioral research_ , 1(2):245–276, 1966.
  * Chen et al. (2020) Shuxiao Chen, Sifan Liu, and Zongming Ma.  Global and individualized community detection in inhomogeneous multilayer networks.  _arXiv preprint arXiv:2012.00933_ , 2020.
  * Chen and Yang (2021) Xiaohui Chen and Yun Yang.  Cutoff for exact recovery of gaussian mixture models.  _IEEE Transactions on Information Theory_ , 67(6):4223–4238, 2021.
  * Chernoff (1952) Herman Chernoff.  A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observations.  _The Annals of Mathematical Statistics_ , pages 493–507, 1952.
  * Dasgupta (2008) Sanjoy Dasgupta.  _The hardness of k-means clustering_.  Department of Computer Science and Engineering, University of California …, 2008.
  * De Lathauwer et al. (2000) Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle.  A multilinear singular value decomposition.  _SIAM journal on Matrix Analysis and Applications_ , 21(4):1253–1278, 2000.
  * Ding et al. (2019) Yunzi Ding, Dmitriy Kunisky, Alexander S Wein, and Afonso S Bandeira.  Subexponential-time algorithms for sparse pca.  _arXiv preprint arXiv:1907.11635_ , 2019.
  * Dong et al. (2012) Xiaowen Dong, Pascal Frossard, Pierre Vandergheynst, and Nikolai Nefedov.  Clustering with multi-layer graphs: A spectral perspective.  _IEEE Transactions on Signal Processing_ , 60(11):5820–5831, 2012.
  * Fei and Chen (2018) Yingjie Fei and Yudong Chen.  Hidden integrality of sdp relaxations for sub-gaussian mixture models.  In _Conference On Learning Theory_ , pages 1931–1965. PMLR, 2018\.
  * Gao and Zhang (2022) Chao Gao and Anderson Y Zhang.  Iterative algorithm for discrete structure recovery.  _The Annals of Statistics_ , 50(2):1066–1094, 2022.
  * Gao et al. (2018) Chao Gao, Zongming Ma, Anderson Y Zhang, and Harrison H Zhou.  Community detection in degree-corrected block models.  _The Annals of Statistics_ , 46(5):2153–2185, 2018.
  * Gao et al. (2021) Xu Gao, Weining Shen, Liwen Zhang, Jianhua Hu, Norbert J Fortin, Ron D Frostig, and Hernando Ombao.  Regularized matrix data clustering and its application to image analysis.  _Biometrics_ , 77(3):890–902, 2021.
  * Gavish and Donoho (2017) Matan Gavish and David L Donoho.  Optimal shrinkage of singular values.  _IEEE Transactions on Information Theory_ , 63(4):2137–2152, 2017.
  * Guo et al. (2010) Jian Guo, Elizaveta Levina, George Michailidis, and Ji Zhu.  Pairwise variable selection for high-dimensional model-based clustering.  _Biometrics_ , 66(3):793–804, 2010.
  * Hajek et al. (2016) Bruce Hajek, Yihong Wu, and Jiaming Xu.  Achieving exact cluster recovery threshold via semidefinite programming.  _IEEE Transactions on Information Theory_ , 62(5):2788–2797, 2016.
  * Han et al. (2015) Qiuyi Han, Kevin Xu, and Edoardo Airoldi.  Consistent estimation of dynamic and multi-layer block models.  In _International Conference on Machine Learning_ , pages 1511–1520. PMLR, 2015.
  * Han et al. (2022a) Rungang Han, Yuetian Luo, Miaoyan Wang, and Anru R Zhang.  Exact clustering in tensor block model: Statistical optimality and computational limit.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 84(5):1666–1698, 2022a.
  * Han et al. (2022b) Rungang Han, Rebecca Willett, and Anru R Zhang.  An optimal statistical and computational framework for generalized tensor estimation.  _The Annals of Statistics_ , 50(1):1–29, 2022b.
  * Hoff et al. (2002) Peter D Hoff, Adrian E Raftery, and Mark S Handcock.  Latent space approaches to social network analysis.  _Journal of the american Statistical association_ , 97(460):1090–1098, 2002.
  * Holland et al. (1983) Paul W Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt.  Stochastic blockmodels: First steps.  _Social networks_ , 5(2):109–137, 1983.
  * Hopkins (2018) Samuel Hopkins.  _Statistical inference and the sum of squares method_.  PhD thesis, Cornell University, 2018.
  * Hu et al. (2020) Wei Hu, Weining Shen, Hua Zhou, and Dehan Kong.  Matrix linear discriminant analysis.  _Technometrics_ , 62(2):196–205, 2020.
  * Huang et al. (2022) Hsin-Hsiung Huang, Feng Yu, Xing Fan, and Teng Zhang.  Robust regularized low-rank matrix models for regression and classification.  _arXiv preprint arXiv:2205.07106_ , 2022.
  * Jin et al. (2016) Chi Jin, Yuchen Zhang, Sivaraman Balakrishnan, Martin J Wainwright, and Michael I Jordan.  Local maxima in the likelihood of gaussian mixture models: Structural results and algorithmic consequences.  _Advances in neural information processing systems_ , 29, 2016.
  * Jin et al. (2017) Jiashun Jin, Zheng Tracy Ke, and Wanjie Wang.  Phase transitions for high dimensional clustering and related problems.  _The Annals of Statistics_ , 45(5):2151–2189, 2017.
  * Jing et al. (2021) Bing-Yi Jing, Ting Li, Zhongyuan Lyu, and Dong Xia.  Community detection on mixture multilayer networks via regularized tensor decomposition.  _The Annals of Statistics_ , 49(6):3181–3205, 2021.
  * Kolda and Bader (2009) Tamara G Kolda and Brett W Bader.  Tensor decompositions and applications.  _SIAM review_ , 51(3):455–500, 2009.
  * Koltchinskii and Lounici (2017) Vladimir Koltchinskii and Karim Lounici.  Concentration inequalities and moment bounds for sample covariance operators.  _Bernoulli_ , 23(1):110–133, 2017.
  * Kumar et al. (2011) Abhishek Kumar, Piyush Rai, and Hal Daume.  Co-regularized multi-view spectral clustering.  _Advances in neural information processing systems_ , 24, 2011.
  * Kumar et al. (2004) Amit Kumar, Yogish Sabharwal, and Sandeep Sen.  A simple linear time (1+/spl epsiv/)-approximation algorithm for k-means clustering in any dimensions.  In _45th Annual IEEE Symposium on Foundations of Computer Science_ , pages 454–462. IEEE, 2004.
  * Kunisky et al. (2019) Dmitriy Kunisky, Alexander S Wein, and Afonso S Bandeira.  Notes on computational hardness of hypothesis testing: Predictions using the low-degree likelihood ratio.  _arXiv preprint arXiv:1907.11636_ , 2019.
  * Larremore et al. (2013) Daniel B Larremore, Aaron Clauset, and Caroline O Buckee.  A network approach to analyzing highly recombinant malaria parasite genes.  _PLoS computational biology_ , 9(10):e1003268, 2013.
  * Levin et al. (2019) Keith Levin, Asad Lodhia, and Elizaveta Levina.  Recovering low-rank structure from multiple networks with unknown edge distributions.  _arXiv preprint arXiv:1906.07265_ , 2019.
  * Li et al. (2010) Bing Li, Min Kyung Kim, and Naomi Altman.  On dimension folding of matrix-or array-valued statistical objects.  _The Annals of Statistics_ , 38(2):1094–1121, 2010.
  * Liu et al. (2022) Tianqi Liu, Ming Yuan, and Hongyu Zhao.  Characterizing spatiotemporal transcriptome of the human brain via low-rank tensor decomposition.  _Statistics in Biosciences_ , pages 1–29, 2022.
  * Lloyd (1982) Stuart Lloyd.  Least squares quantization in pcm.  _IEEE transactions on information theory_ , 28(2):129–137, 1982.
  * Löffler et al. (2020) Matthias Löffler, Alexander S Wein, and Afonso S Bandeira.  Computationally efficient sparse clustering.  _arXiv preprint arXiv:2005.10817_ , 2020.
  * Löffler et al. (2021) Matthias Löffler, Anderson Y Zhang, and Harrison H Zhou.  Optimality of spectral clustering in the gaussian mixture model.  _The Annals of Statistics_ , 49(5):2506–2530, 2021.
  * Lu and Zhou (2016) Yu Lu and Harrison H Zhou.  Statistical and computational guarantees of lloyd’s algorithm and its variants.  _arXiv preprint arXiv:1612.02099_ , 2016.
  * Luo and Zhang (2022) Yuetian Luo and Anru R Zhang.  Tensor clustering with planted structures: Statistical optimality and computational limits.  _The Annals of Statistics_ , 50(1):584–613, 2022\.
  * Lyu and Xia (2022) Zhongyuan Lyu and Dong Xia.  Optimal estimation and computational limit of low-rank gaussian mixtures.  _arXiv preprint arXiv:2201.09040_ , 2022.
  * Lyu et al. (2021) Zhongyuan Lyu, Dong Xia, and Yuan Zhang.  Latent space model for higher-order networks and generalized tensor decomposition.  _arXiv preprint arXiv:2106.16042_ , 2021.
  * Mahajan et al. (2009) Meena Mahajan, Prajakta Nimbhorkar, and Kasturi Varadarajan.  The planar k-means problem is np-hard.  In _International workshop on algorithms and computation_ , pages 274–285. Springer, 2009.
  * Mai et al. (2021) Qing Mai, Xin Zhang, Yuqing Pan, and Kai Deng.  A doubly enhanced em algorithm for model-based tensor clustering.  _Journal of the American Statistical Association_ , pages 1–15, 2021\.
  * Mendelson (2016) Shahar Mendelson.  Upper bounds on product and multiplier empirical processes.  _Stochastic Processes and their Applications_ , 126(12):3652–3680, 2016.
  * Ndaoud (2018) Mohamed Ndaoud.  Sharp optimal recovery in the two-component gaussian mixture model.  _arXiv preprint arXiv:1812.08078_ , 2018.
  * Richard and Montanari (2014) Emile Richard and Andrea Montanari.  A statistical model for tensor pca.  _Advances in neural information processing systems_ , 27, 2014.
  * Stanley et al. (2016) Natalie Stanley, Saray Shai, Dane Taylor, and Peter J Mucha.  Clustering network layers with the strata multilayer stochastic block model.  _IEEE transactions on network science and engineering_ , 3(2):95–105, 2016.
  * Sun and Li (2019) Will Wei Sun and Lexin Li.  Dynamic tensor clustering.  _Journal of the American Statistical Association_ , 114(528):1894–1907, 2019.
  * Vempala and Wang (2004) Santosh Vempala and Grant Wang.  A spectral algorithm for learning mixture models.  _Journal of Computer and System Sciences_ , 68(4):841–860, 2004.
  * Verzelen and Arias-Castro (2017) Nicolas Verzelen and Ery Arias-Castro.  Detection and feature selection in sparse mixture models.  _The Annals of Statistics_ , 45(5):1920–1950, 2017.
  * Wang and Zeng (2019) Miaoyan Wang and Yuchen Zeng.  Multiway clustering via tensor block models.  _Advances in neural information processing systems_ , 32, 2019.
  * Witten and Tibshirani (2010) Daniela M Witten and Robert Tibshirani.  A framework for feature selection in clustering.  _Journal of the American Statistical Association_ , 105(490):713–726, 2010.
  * Wu and Zhou (2019) Yihong Wu and Harrison H Zhou.  Randomly initialized em algorithm for two-component gaussian mixture achieves near optimality in $o(\sqrt{n})$ iterations.  _arXiv preprint arXiv:1908.10935_ , 2019.
  * Xia (2021) Dong Xia.  Normal approximation and confidence region of singular subspaces.  _Electronic Journal of Statistics_ , 15(2):3798–3851, 2021.
  * Xia and Yuan (2021) Dong Xia and Ming Yuan.  Statistical inferences of linear forms for noisy matrix completion.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 83(1):58–77, 2021.
  * Xia and Zhou (2019) Dong Xia and Fan Zhou.  The sup-norm perturbation of hosvd and low rank tensor denoising.  _The Journal of Machine Learning Research_ , 20(1):2206–2247, 2019.
  * Xia et al. (2022) Dong Xia, Anru R Zhang, and Yuchen Zhou.  Inference for low-rank tensors—no need to debias.  _The Annals of Statistics_ , 50(2):1220–1245, 2022.
  * Zhang and Zhou (2022) Anderson Y Zhang and Harrison H Zhou.  Leave-one-out singular subspace perturbation analysis for spectral clustering.  _arXiv preprint arXiv:2205.14855_ , 2022.
  * Zhang and Xia (2018) Anru Zhang and Dong Xia.  Tensor svd: Statistical and computational limits.  _IEEE Transactions on Information Theory_ , 64(11):7311–7338, 2018.
  * Zhang et al. (1995) Xiao Lei Zhang, Henri Begleiter, Bernice Porjesz, Wenyu Wang, and Ann Litke.  Event related potentials during object recognition tasks.  _Brain research bulletin_ , 38(6):531–538, 1995\.
  * Zhou and Li (2014) Hua Zhou and Lexin Li.  Regularized matrix regression.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(2):463–483, 2014.

##  Appendix A Proofs of Main Theorems

Throughout the proofs, we use $c,C,C^{\prime}$ to represent generic absolute
constants, whose actual values may vary in different formulas.

###  A.1 Proof of Theorem 1

##### Step 1: Notations and Good Initialization

We need to introduce some notations to simplify the presentation of our proof.
Recall the individual signal strength is defined as

$\lambda=\min_{k\in[K]}\sigma_{\textsf{\tiny min}}(\mathbf{M}_{k})$

Note in our setting, we simply have
$\lambda\gtrsim\kappa_{0}^{-1}r^{-1/2}\max_{a\neq
b}\left\|\mathbf{M}_{a}-\mathbf{M}_{b}\right\|_{\rm
F}\geq\kappa_{0}^{-1}r^{-1/2}\Delta$.
Define the frobenius error with respect to the true label $\mathbf{s}^{\ast}$:

$\ell(\mathbf{s},\mathbf{s}^{\ast}):=\sum_{i=1}^{n}\left\|\mathbf{M}_{s_{i}}-\mathbf{M}_{s_{i}^{\ast}}\right\|_{\rm
F}^{2}$

as well as the corresponding hamming loss:

$h(\mathbf{s},\mathbf{s}^{\ast}):=\sum_{i=1}^{n}\mathbb{I}\left(s_{i}\neq
s_{i}^{\ast}\right)$

A simple relation is that
$h(\mathbf{s},\mathbf{s}^{\ast})\leq\Delta^{-2}\cdot\ell(\mathbf{s},\mathbf{s}^{\ast})$
due to the fact

$\sum_{i=1}^{n}\|\mathbf{M}_{s_{i}}-\mathbf{M}_{s_{i}^{\ast}}\|_{\rm
F}^{2}\geq\sum_{i=1}^{n}\mathbb{I}\left(s_{i}\neq
s_{i}^{\ast}\right)\Delta^{2}.$

Note that, by definition
$\ell_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})=\sum_{i=1}^{n}\left\|\mathbf{M}_{s_{i}^{(0)}}-\mathbf{M}_{\pi(s_{i}^{\ast})}\right\|_{\rm
F}^{2}$ for some permutation $\pi$, we can always relabel our
$\mathbf{M}_{1},\cdots,\mathbf{M}_{K}$ to
$\mathbf{M}_{\pi(1)},\cdots,\mathbf{M}_{\pi(K)}$ after initialization.
Therefore, without loss of generality we can assume $\pi=\text{Id}$ and hence
$\ell(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})=\ell_{\textsf{c}}(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})$.
As a result of condition (10), we also have

$\displaystyle
h(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})\leq\frac{\ell(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})}{\Delta^{2}}=o\left(\frac{\alpha
n}{\kappa_{0}^{2}K}\right)$ (20)

Note that (10) can be equivalently expressed as
$\ell(\widehat{}\mathbf{s}^{(0)},\mathbf{s}^{\ast})\leq\tau$ for some
$\tau=o\left(\kappa_{0}^{-2}{\alpha n\Delta^{2}}/{K}\right)$ and hence
$\Delta^{2}\gg\kappa_{0}^{2}K\tau/(\alpha n)$.

##### Step 2: Iterative Convergence

We then analyze the convergence property of low-rank Lloyd algorithm. Without
loss of generality, given the labelling $\widehat{}\mathbf{s}^{(t-1)}$ at the
$(t-1)$-th iteration, we investigate the behavior of
$\widehat{}\mathbf{s}^{(t)}$, i.e., after one iteration of Lloyd algorithm.