# Optimal Variable Clustering for High-Dimensional Matrix Valued Data

Inbeom Lee Siyi Deng Yang Ning  Department of Statistics and Data Science,
Cornell University, Ithaca, NY. E-mail<EMAIL_ADDRESS>of
Statistics and Data Science, Cornell University, Ithaca, NY. E-mail:
<EMAIL_ADDRESS>of Statistics and Data Science, Cornell
University, Ithaca, NY. E-mail<EMAIL_ADDRESS>
###### Abstract

Matrix valued data has become increasingly prevalent in many applications.
Most of the existing clustering methods for this type of data are tailored to
the mean model and do not account for the dependence structure of the
features, which can be very informative, especially in high-dimensional
settings. To extract the information from the dependence structure for
clustering, we propose a new latent variable model for the features arranged
in matrix form, with some unknown membership matrices representing the
clusters for the rows and columns. Under this model, we further propose a
class of hierarchical clustering algorithms using the difference of a weighted
covariance matrix as the dissimilarity measure. Theoretically, we show that
under mild conditions, our algorithm attains clustering consistency in the
high-dimensional setting. While this consistency result holds for our
algorithm with a broad class of weighted covariance matrices, the conditions
for this result depend on the choice of the weight. To investigate how the
weight affects the theoretical performance of our algorithm, we establish the
minimax lower bound for clustering under our latent variable model. Given
these results, we identify the optimal weight in the sense that using this
weight guarantees our algorithm to be minimax rate-optimal in terms of the
magnitude of some cluster separation metric. The practical implementation of
our algorithm with the optimal weight is also discussed. Finally, we conduct
simulation studies to evaluate the finite sample performance of our algorithm
and apply the method to a genomic dataset.

Keywords: Clustering, matrix data, high dimensional estimation, minimax
optimality, latent variable model, hierarchical algorithm

##  1 Introduction

Cluster analysis is one of the most important unsupervised learning techniques
and has been widely used to discover the underlying group structure in data,
arising in many applications including economics, image analysis, psychology
and the biomedical sciences (Everitt et al., 2011; Kogan, 2007). In these
applications, matrix valued data is becoming increasingly prevalent. For
example, in genetic studies, one may observe data matrices of dimension
$p\times q$ from $n$ subjects, where the $(j,k)$th entry corresponds to the
expression value of the $j$th gene at the $k$th tissue (Zahn et al., 2007).
The biologist is often interested in identifying clusters of genes that share
similar biological functions and also clusters of similar tissues. Similarly,
in functional magnetic resonance imaging (fMRI) studies, to understand how the
brain connectivity structure changes under different tasks/stimuli,
researchers can measure the blood oxygen level (BOLD) within each region of
interest (ROI) from the brain under a variety of conditions (Mitchell et al.,
2008). The data points from each participant can be stacked as a matrix, in
which the rows and columns correspond to different ROIs and tasks/stimuli,
respectively. Given the data from $n$ participants, it is of interest to
simultaneously cluster the ROIs and tasks/stimuli. Such clustering results can
be used as a dimension reduction step to further investigate brain
connectivity networks (Eisenach et al., 2020). Driven by these applications,
the goal of this paper is to develop a statistical clustering framework for
matrix valued data with theoretical guarantees. In particular, assuming that
$n$ i.i.d. samples $X^{(1)},...,X^{(n)}$ of a random matrix
$X\in\mathbb{R}^{p\times q}$ are observed, we would like to recover the
cluster membership of the rows and the columns of the feature matrix $X$.

In the literature, clustering a data matrix is often known as biclustering
(Hartigan, 1972; Madeira and Oliveira, 2004). Most of the existing
biclustering methods can be classified into the following categories: (1)
hierarchical approaches based on the dendrogram (Hastie et al., 2009); (2)
extensions of K-means (Fraiman and Li, 2022); (3) the penalized likelihood
approach (Tan and Witten, 2014); (4) convex clustering via fused lasso (Chi et
al., 2017); and (5) clustering based on the singular value decomposition (Sill
et al., 2011; Lee et al., 2010). In these existing works, the goal is to
simultaneously cluster $n$ samples and $d$ features stacked as an $n\times d$
data matrix, whereas in our problem we are interested in clustering a $p\times
q$ feature matrix. Since the features are often correlated, the $p\times q$
feature matrix induces a sophisticated but informative dependence structure
that is not accounted for in the existing biclustering methods. Thus, these
existing biclustering methods may not be appropriate in our setting. More
recently, multiway clustering, also known as tensor clustering, is attracting
increasing attention (Mankad and Michailidis, 2014; Zhao et al., 2016; Chi et
al., 2020; Sun and Li, 2019; Wang et al., 2019; Wang and Zeng, 2019). Since we
observe $n$ i.i.d $p\times q$ feature matrices, we can view our data as a
3-way tensor with dimension $p\times q\times n$. The existing tensor
clustering methods such as the tensor block model (Wang and Zeng, 2019; Chi et
al., 2020) aim to cluster each mode of the tensor from the mean model. While
these approaches enjoy great success in many applications and can be adapted
to our setting, they may not perform well when the dependence structure of the
features is informative. There are several recent works on tensor clustering
that utilize the dependence structure (Deng and Zhang (2022), Mai et al.
(2022)). However, these methods focus on clustering over the observations and
not over the features (i.e., the rows and the columns of $X$).

In this paper, our first contribution is to propose a new latent variable
model for clustering matrix valued data. Assume that the features are stacked
as a random matrix $X\in\mathbb{R}^{p\times q}$, which follows

$X=AZB^{T}+\Gamma,$ (1.1)

where without loss of generality, $Z\in\mathbb{R}^{K_{1}\times K_{2}}$ is a
latent variable matrix with $\mathbb{E}(Z)=0$. $A\in\mathbb{R}^{p\times
K_{1}}$ and $B\in\mathbb{R}^{q\times K_{2}}$ are the unknown binary membership
matrices for the rows and columns, respectively, and
$\Gamma\in\mathbb{R}^{p\times q}$ represents the random noise matrix with
entries that have mean 0 and variance $\sigma_{ab}^{2}=\operatorname{{\rm
Var}}(\Gamma_{ab})$. We assume that the entries in the noise matrix $\Gamma$
are mutually independent and are also independent of the entries in $Z$.
Entries of the membership matrix $A$ take values in $\\{0,1\\}$, such that
$A_{ak}=1$ if row $a$ belongs to row cluster $k$ and $A_{ak}=0$ otherwise. In
this paper, we focus on the non-overlapping and exhaustive clustering
scenario. That is, for each row $a$, there exists one and only one cluster $k$
with $A_{ak}=1$. The same requirement holds for the membership matrix $B$. To
see why $A$ and $B$ are interpreted as membership matrices, we note that, for
any feature $X_{ab}$, if $A_{ak}=1$ and $B_{b\ell}=1$ for some $k$ and $\ell$,
then model (1.1) implies $X_{ab}=Z_{k\ell}+\Gamma_{ab}$. That is, it implies
that the feature $X_{ab}$ is associated with the latent variable $Z_{k\ell}$.
For this reason, we say that $X_{ab}$ belongs to the row cluster $k$ and
column cluster $\ell$. Under model (1.1), we can formally define the row
clusters as the partition

$\mathcal{G}^{(r)}=\\{G_{1}^{(r)},...,G_{K_{1}}^{(r)}\\}\leavevmode\nobreak\
\leavevmode\nobreak\ \textrm{of}\leavevmode\nobreak\ \leavevmode\nobreak\
\\{1,...,p\\},\leavevmode\nobreak\ \leavevmode\nobreak\
\textrm{where}\leavevmode\nobreak\ \leavevmode\nobreak\
G_{k}^{(r)}=\\{a:A_{ak}=1\\}$ (1.2)

for any $1\leq k\leq K_{1}$, and $K_{1}$ is the unknown number of row
clusters. When the context is clear, $\mathcal{G}$ will be used for notational
simplicity. Without loss of generality, we focus on how to recover the unknown
membership matrix $A$ for the rows (or equivalently $\mathcal{G}^{(r)}$) up to
label switching.

Model (1.1) can be viewed as the extension of the G-block model for clustering
a random vector in Bunea et al. (2020) to matrix valued data. Indeed, if we
vectorize the matrix $X$, model (1.1) is equivalent to
$\textrm{vec}(X)=M\textrm{vec}(Z)+\textrm{vec}(\Gamma)$ with $M=B\otimes A$,
where $\textrm{vec}(X)$ denotes the vectorization of $X$, formed by stacking
the columns of $X$ into a single column vector, and $\otimes$ denotes the
Kronecker product. Thus, compared to Bunea et al. (2020) which allows $M$ to
be any unstructured $(pq)\times(K_{1}K_{2})$ membership matrix, we impose the
Kronecker product structure to the membership matrix $M$. While our model for
$\textrm{vec}(X)$ is more restrictive than the model in Bunea et al. (2020),
it actually comes with two advantages for matrix clustering. First, as seen
above, $A$ and $B$ are interpreted as the membership matrices for the rows and
columns. Ignoring the Kronecker product structure and directly applying the
model in Bunea et al. (2020) would no longer produce interpretable results for
matrix clustering, as shown in Figure 3 in Section A in the Supplementary
Material. Second, the Kronecker product of $A$ and $B$ provides a more
parsimonious parametrization for the unknown membership matrix $M$, leading to
stronger theoretical guarantees on clustering.

It is also worth mentioning that our proposed latent variable model shares
similarities with the stochastic block model (SBM) widely used in community
detection in network analysis. First introduced in Holland et al. (1983), the
stochastic block model assumes that nodes of a network are partitioned into
subgroups called blocks and the distribution of the ties between nodes is
dependent on the blocks to which the nodes belong. More specifically, they
imposed the model $\mathbb{E}(A)=ZBZ^{T}$ where $A$ is the $n\times n$
adjacency matrix, $Z$ is the $n\times K$ binary block membership matrix and
$B$ is the $K\times K$ matrix whose element $B_{ij}$ represents the
probability that a node from block $i$ is connected to a node in block $j$. We
defer the details on further development of the method to Airoldi et al.
(2008); Zhang et al. (2020); Abbe (2017) and the references therein. It must
be noted, however, that community detection with stochastic block type models
is inherently different from our setting since the former is focused on
clustering the nodes (i.e., the samples), whereas our matrix clustering
setting is focused on clustering the features (the rows and the columns of
$X$).

To recover the unknown membership matrix $A$ (and similarly for $B$), our
second contribution is to propose a class of hierarchical clustering methods
based on the weighted covariance matrix $\Sigma_{p,W}=\mathbb{E}(XWX^{T})$ for
some positive semi-definite matrix $W\in\mathbb{R}^{q\times q}$ to be chosen.
We use the difference of this weighted covariance matrix as the dissimilarity
measure in our hierarchical algorithm to recover the membership matrix $A$. To
establish the theoretical guarantees of our algorithm, we introduce the metric
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$ to quantify how well the clusters are
separated. The precise definition is detailed in Section 3.1. Theoretically,
we develop a general result on the clustering consistency of our hierarchical
algorithm that holds for a broad class of weight matrices $W$. To attain
clustering consistency for the rows, we require the cluster separation metric
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$ to be no smaller than the order of
$(\frac{\log p}{nK_{2}})^{1/2}$, where $K_{2}$ is the unknown number of column
clusters. The implication is that with the help of a larger $K_{2}$,
clustering the rows of $X$ becomes easier.

To investigate the optimality of our algorithm, we establish the minimax lower
bound for clustering under our latent variable model. While the clustering
consistency property holds for a broad class of weight matrices $W$, our
algorithm with a generic weight $W$ may not be minimax optimal. To derive an
optimal clustering algorithm for the rows, the key is to account for the
information in the column clusters. The intuition is that with a more accurate
column clustering result (i.e., estimate of $B$), we can decorrelate the
dependence structure of $X$ and reduce the noise induced by $\Gamma$, which in
turn improves the clustering accuracy for the rows. Following this argument,
we define the optimal weight $W_{O}$, and propose to estimate it by
$\widehat{W}_{O}=\widehat{B}(\widehat{B}^{T}\widehat{B})^{-2}\widehat{B}^{T}/s$,
where $\widehat{B}$ is an estimate of the column membership matrix $B$ and $s$
denotes the estimated number of clusters from $\widehat{B}$. Under mild
conditions, we show that the proposed algorithm with the estimated optimal
weight attains the minimax lower bound, and therefore is rate-optimal for
clustering. To the best of our knowledge, our paper is the first to formally
establish minimax optimality for clustering matrix valued data. From a
technical perspective, to show the optimality of our algorithm, the main
challenge is to quantify the “stability” of the algorithm with respect to an
imperfect estimate $\widehat{B}$ in the weight matrix. We provide sufficient
conditions to show the stability of the algorithm under some additional
modeling assumptions. In practice, the hierarchical algorithm can be applied
iteratively to cluster the rows and columns. Finally, we conduct extensive
numerical studies to support our theoretical results.

The rest of the paper is organized as follows. In Section 2, we propose the
hierarchical algorithm and discuss the advantage of using the optimal weight,
$W_{O}$. In Section 3, we define a notion of cluster separation and use it to
establish the clustering consistency of our algorithm and to derive the
minimax lower bound. Since the algorithm using the optimal weight depends on
the initial estimate $\widehat{B}$, in Section 4, we verify the cluster
separation condition and stability condition required in Theorem 3.1 for
clustering consistency in matrix normal models. The practical implementation
of the algorithm is discussed in Section 5. Three further extensions of our
method - the dependent noise model, the nested clustering method to
incorporate both mean and covariance information, and higher order tensor
models - are discussed in Section 6. The simulation results and real data
analysis are presented in Sections 7 and 8, respectively. The paper concludes
with a discussion in Section 9.

Notation. For any $1\leq a\neq b\leq p$, we write
$a\underset{\mathcal{G}}{\sim}b$ if $a$ and $b$ are in the same cluster (i.e.,
$a,b\in G_{k}$ for some $k$). Otherwise, we write
$a\underset{\mathcal{G}}{\nsim}b$. For a matrix $X$, we use the following
norms: $||X||_{\mathrm{max}}=\max_{i,j}|X_{ij}|$,
$||X||_{\infty}=\max_{i}\sum_{j=1}^{q}|X_{ij}|$,
$||X||_{F}=\sqrt{\sum_{i,j}X_{ij}^{2}}$. $||X||_{\text{op}}$ denotes the
largest singular value of $X$. The largest and smallest eigenvalues are
denoted by $\lambda_{\max}(\cdot)$ and $\lambda_{\min}(\cdot)$. We use
$X_{\cdot j}$ and $X_{j\cdot}$ to denote the $j$th column and row of $X$,
respectively. For two positive sequences $a_{n}$ and $b_{n}$, we write
$a_{n}\asymp b_{n}$ if $C\leq a_{n}/b_{n}\leq C^{\prime}$ for all $n$ for some
constant $C,C^{\prime}>0$. Similarly, we use $a_{n}\lesssim b_{n}$
($a_{n}\gtrsim b_{n}$) to denote $a_{n}\leq Cb_{n}$ ($a_{n}\geq Cb_{n}$) for
all $n$ for some constant $C>0$.

##  2 Methodology

###  2.1 Hierarchical Clustering via Weighted Covariance Differences

Recall that the random matrix $X$ follows model (1.1). In this section, we
propose a class of clustering methods for the rows of $X$ based on the
weighted covariance matrix $\Sigma_{p,W}=\mathbb{E}(XWX^{T})$, where
$W\in\mathbb{R}^{q\times q}$ is some positive semi-definite matrix to be
chosen. We add a subscript $p$ to indicate $\Sigma_{p,W}$ is a $p\times p$
matrix corresponding to the rows. The same type of method can be used to
cluster the columns of $X$.

In the following, we first outline how to identify the unknown membership
matrix $A$ from the weighted covariance matrix $\Sigma_{p,W}$ on the
population level. Under model (1.1), by the independence between $Z$ and
$\Gamma$, we obtain that

$\Sigma_{p,W}=A\mathbb{E}(ZB^{T}WBZ^{T})A^{T}+\mathbb{E}(\Gamma W\Gamma^{T}).$
(2.1)

To recover the membership matrix $A$ from $\Sigma_{p,W}$, one needs to first
separate the two matrices $A\mathbb{E}(ZB^{T}WBZ^{T})A^{T}$ and
$\mathbb{E}(\Gamma W\Gamma^{T})$. Noting that $\mathbb{E}(\Gamma W\Gamma^{T})$
is a diagonal matrix as the elements in $\Gamma$ are mutually independent and
the weight matrix $W$ is deterministic, we therefore focus on the non-diagonal
entries of $\Sigma_{p,W}$, that is
$[\Sigma_{p,W}]_{ac}=[A\mathbb{E}(ZB^{T}WBZ^{T})A^{T}]_{ac}$ for any $1\leq
a\neq c\leq p$. By the definition of the membership matrix $A$, for any $a\in
G_{k}$, $c\in G_{\ell}$ and $a\neq c$, we have

$\big{(}\Sigma_{p,W}\big{)}_{ac}=\Big{[}A\mathbb{E}(ZB^{T}WBZ^{T})A^{T}\Big{]}_{ac}=\Big{[}\mathbb{E}(ZB^{T}WBZ^{T})\Big{]}_{k\ell}.$
(2.2)

In view of (2.2), the within-cluster covariance difference
$(\Sigma_{p,W})_{ac}-(\Sigma_{p,W})_{bc}$ with
$a\underset{\mathcal{G}}{\sim}b$ is always 0 for any $c\neq a,b$. In other
words, as long as $(\Sigma_{p,W})_{ac}-(\Sigma_{p,W})_{bc}$ is nonzero for
some $c$, it indicates that $a$ and $b$ are not in the same cluster. Thus, the
covariance difference is indicative of the clustering structure. Following
Bunea et al. (2020), we formally define the covariance difference (COD) as

$\displaystyle\mathrm{COD}_{\Sigma_{p,W}}(a,b)$
$\displaystyle:=\underset{c\neq a,b}{\text{max}}\leavevmode\nobreak\
\Big{|}\big{(}\Sigma_{p,W}\big{)}_{ac}-\big{(}\Sigma_{p,W}\big{)}_{bc}\Big{|}.$
(2.3)

From the above argument, we have $\mathrm{COD}_{\Sigma_{p,W}}(a,b)=0$ if
$a\underset{\mathcal{G}}{\sim}b$. Moreover, if
$\mathrm{COD}_{\Sigma_{p,W}}(a,b)>0$ holds for all
$a\underset{\mathcal{G}}{\nsim}b$, we are able to identify all the clusters.
Let

$\displaystyle\mathrm{MCOD}(\Sigma_{p,W})$
$\displaystyle:=\min_{a\underset{\mathcal{G}}{\nsim}b}\leavevmode\nobreak\
\text{COD}_{\Sigma_{p,W}}(a,b)$ (2.4)

denote the minimum COD value over all possible
$a\underset{\mathcal{G}}{\nsim}b$. On the population level, provided
$\mathrm{MCOD}(\Sigma_{p,W})>0$, the membership matrix $A$ (or equivalently
$\mathcal{G}^{(r)}$ in (1.2)) is identifiable from the weighted covariance
matrix $\Sigma_{p,W}$ up to label switching.

Based on the above results on the population level, we will develop a
hierarchical clustering algorithm to estimate the membership matrix $A$ (or
equivalently $\mathcal{G}^{(r)}$ in (1.2)). Given $n$ i.i.d. samples
$X^{(1)},...,X^{(n)}$, we first estimate $\Sigma_{p,W}$ with

$\widehat{\Sigma}_{p,W}=\frac{1}{n}\sum_{i=1}^{n}X^{(i)}WX^{(i)T},$

and then plug this into (2.3) to form a dissimilarity measure
$\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(a,b)$ for any $1\leq a,b\leq p$. The
hierarchical algorithm starts with every variable representing a singleton
cluster. At each step, the closest two clusters are merged into one single
cluster based on the following dissimilarity measure between two sets $I$ and
$J$

$\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(I,J)=\max_{a\in I,b\in
J}\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(a,b).$ (2.5)

We refer to Hastie et al. (2009) for alternative definitions of dissimilarity
between two clusters and further discussions. Finally, we terminate this
process and report the clusters when the dissimilarity measure
$\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(I,J)$ exceeds a threshold value
$\alpha>0$. This hierarchical algorithm, summarized in Algorithm 1, improves
the existing COD algorithm proposed by Bunea et al. (2020). First, our
algorithm satisfies the so-called monotonicity property, which states that
merged clusters always have smaller values of
$\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(I,J)$ than unmerged ones (Hastie et
al., 2009). However, the algorithm in Bunea et al. (2020) may merge indices
$a$ and $b$ into one cluster even if there exists another index $c$ with
$\mathrm{COD}_{\Sigma_{p,W}}(a,c)$ smaller than
$\mathrm{COD}_{\Sigma_{p,W}}(a,b)$. We include in Section B in the
Supplementary Material a toy example that illustrates this point. Empirically,
we find that our hierarchical algorithm produces more stable clustering
results than the algorithm in Bunea et al. (2020). Second, the hierarchical
algorithm is more flexible in incorporating side information, such as the
number of row clusters $K_{1}$. While our Algorithm 1 does not require the
user to know $K_{1}$, with such information from domain knowledge or existing
literature, the algorithm is expected to yield more reliable clustering
results.

INPUT: The estimated weighted covariance matrix $\widehat{\Sigma}_{p,W}$ and a
threshold $\alpha>0$.

  * (1)

Calculate $\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(a,b)$ for $1\leq a,b\leq p$.

  * (2)

Create a hierarchical tree based on the value of

    * $\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(I,K)=\underset{a\in I,\leavevmode\nobreak\ b\in K}{\max}\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(a,b)$ for sets $I$ and $K$.

  * (3)

Use the threshold value $\alpha$ to cut the tree and obtain the estimated row
clustering $\widehat{\mathcal{G}}^{(r)}$. More precisely, we use the following
rule to find the clusters:

    * For any two sets of candidate clusters $I$ and $K$ from the hierarchical tree,

merge them into one cluster if and only if
$\mathrm{COD}_{\widehat{\Sigma}_{p,W}}(I,K)\leq\alpha$.

Algorithm 1 Hierarchical Algorithm with COD

###  2.2 Optimal Choice of $W$

While our Algorithm 1 can be applied with any weight matrix $W$ in
$\widehat{\Sigma}_{p,W}$, the empirical and theoretical performance of the
algorithm critically depends on the choice of $W$. In practice, the simplest
choice of $W$ could be $W_{I}=I_{q}/q$, where $I_{q}$ is a $q\times q$
identity matrix. With this choice of $W_{I}$, the weighted covariance matrix
$\Sigma_{p,W_{I}}=\frac{1}{q}\sum_{j=1}^{q}\mathbb{E}(X_{\cdot j}X_{\cdot
j}^{T})$ can be interpreted as the average of the second order moment of the
columns of $X$. This “naive” weight $W_{I}$ can be used directly in our
Algorithm 1 (named NAIVE COD) or it can be used as an initial value in a
multi-step iterative algorithm (1-STEP COD, 2-STEP COD), which will be further
discussed in Section 5.

To motivate the development of the optimal choice of $W$, we temporarily
assume that the true column cluster structure (i.e., the membership matrix
$B$) is known up to label switching. We define
$X^{*}=XB(B^{T}B)^{-1}\in\mathbb{R}^{p\times K_{2}}$, which can be interpreted
as the average of $X$ over columns in the same column cluster. To see this,
let us consider a toy example. Assume that $X$ has $q=4$ columns with
$K_{2}=2$ column clusters, where the first two columns belong to cluster 1 and
the last two columns belong to cluster 2. In this case, the membership matrix
$B$ can be written as

$B=\begin{bmatrix}1&0\\\ 1&0\\\ 0&1\\\ 0&1\end{bmatrix}\in\mathbb{R}^{4\times
2}.\leavevmode\nobreak\ \leavevmode\nobreak\ \textrm{Then we
have}\leavevmode\nobreak\ \leavevmode\nobreak\
X^{*}=\begin{bmatrix}\frac{X_{11}+X_{12}}{2}&\frac{X_{13}+X_{14}}{2}\\\ :&:\\\
\frac{X_{p1}+X_{p2}}{2}&\frac{X_{p3}+X_{p4}}{2}\\\
\end{bmatrix}\in\mathbb{R}^{p\times 2}.$

Clearly, the two columns of $X^{*}$ represent the averages of $X$ in the same
column cluster. Inspired by the interpretation of $\Sigma_{p,W_{I}}$, we now
compute the average of the second order moment of the columns of $X^{*}$ as
follows:

$\frac{1}{K_{2}}\sum_{j=1}^{K_{2}}\mathbb{E}(X^{*}_{\cdot j}X_{\cdot
j}^{*T})=\frac{1}{K_{2}}\mathbb{E}(X^{*}X^{*T})=\mathbb{E}(XW_{O}X^{T})=\Sigma_{p,W_{O}},$
(2.6)

where we set $W_{O}=B(B^{T}B)^{-2}B^{T}/K_{2}$ by the definition of $X^{*}$.
This matrix $W_{O}$ is the optimal weight for reasons that will be explained
in Section 3. Intuitively, the weighted covariance matrix $\Sigma_{p,W_{O}}$
with $W_{O}$ is more informative for clustering than $\Sigma_{p,W_{I}}$, as
the random noise in $\Gamma$ is reduced when we construct $X^{*}$ by
aggregating the columns of $X$ in the same column cluster. To better
illustrate this point, we consider a special case. Assume that the columns of
$X$ have $K_{2}$ clusters with equal size $q/K_{2}$. After some algebra, it is
shown that

$\Sigma_{p,W_{I}}=\frac{1}{K_{2}}A\mathbb{E}(ZZ^{T})A^{T}+\frac{1}{q}\mathbb{E}(\Gamma\Gamma^{T}),$

and

$\Sigma_{p,W_{O}}=\frac{1}{K_{2}}A\mathbb{E}(ZZ^{T})A^{T}+\frac{K_{2}}{q}\bigg{\\{}\frac{1}{q}\mathbb{E}(\Gamma\Gamma^{T})\bigg{\\}}.$

Clearly, both $\Sigma_{p,W_{I}}$ and $\Sigma_{p,W_{O}}$ contain the same
amount of row cluster information via the term
$\frac{1}{K_{2}}A\mathbb{E}(ZZ^{T})A^{T}$. However, compared to
$\Sigma_{p,W_{I}}$, the error matrix induced by the covariance of $\Gamma$ is
further reduced by a factor of $K_{2}/q$ in $\Sigma_{p,W_{O}}$. Therefore, we
expect the clustering algorithm using $\Sigma_{p,W_{O}}$ to outperform the
naive method using $\Sigma_{p,W_{I}}$ for recovering the membership matrix
$A$, which is indeed the case in simulations; see Section 7.

Since $W_{O}=B(B^{T}B)^{-2}B^{T}/K_{2}$ depends on the unknown membership
matrix $B$, Algorithm 1 is not directly applicable. In principle, if an
initial estimate of $B$, say $\widehat{B}$, is available, we can plug in the
estimator $\widehat{B}$ and apply Algorithm 1 with
$\widehat{\Sigma}_{p,\widehat{W}_{O}}$, where
$\widehat{W}_{O}=\widehat{B}(\widehat{B}^{T}\widehat{B})^{-2}\widehat{B}^{T}/s$
and $s$ denotes the estimated number of clusters from $\widehat{B}$.
Theoretically, in the next section, we will establish a general result on the
clustering consistency of Algorithm 1 with a data dependent weight matrix
$\widehat{W}$, which covers the case with $\widehat{W}_{O}$. In practice, we
recommend using an iterative hierarchical algorithm to repeatedly cluster the
rows and columns of $X$. The detailed implementation of the algorithm using
the optimal weight is discussed in Section 5.

##  3 Theoretical Guarantees

In this section we establish the theoretical results for our proposed
clustering method. In Section 3.1, we present a general result on the
clustering consistency of Algorithm 1 with a data dependent weight matrix
$\widehat{W}$. Subsequently, we develop the minimax lower bound for the matrix
clustering problem in Section 3.2. In particular, these results imply that
Algorithm 1 with the estimated optimal weight matrix $\widehat{W}_{O}$ defined
in Section 2.2 is minimax optimal for clustering.

###  3.1 Clustering Consistency of Algorithm 1

To formally study the clustering consistency property, we need to first define
a proper notion of cluster separation distance. Recall from the argument in
Section 2.1 that $\mathrm{MCOD}(\Sigma_{p,W})>0$ implies the identifiability
of $A$ up to label switching. One may attempt to use
$\mathrm{MCOD}(\Sigma_{p,W})$ to measure cluster separation. But, using
$\mathrm{MCOD}(\Sigma_{p,W})$ alone is not ideal, as
$\mathrm{MCOD}(\Sigma_{p,W})$ is not invariant to the scale of $W$. To be
specific, we have
$\mathrm{MCOD}(\Sigma_{p,tW})=t\cdot\mathrm{MCOD}(\Sigma_{p,W})$ for any
$t>0$, implying that the MCOD value can be arbitrarily large by rescaling $W$.

In order to resolve this issue, we define $||X||_{W}$ as follows:

$\|X\|_{W}=K_{2}^{1/2}\max_{1\leq a\leq p}\big{\|}L^{T}\operatorname{{\rm
Var}}(X_{a\cdot})L\big{\|}_{F},$

where $L$ is a matrix satisfying $LL^{T}=W$. We note that $\|X\|_{W}$ can be
interpreted as the amount of the variance of the row vector in $X$ reweighted
by $L$. From a technical perspective, such a quantity plays a natural role
when applying the concentration inequality to control
$\widehat{\Sigma}_{p,W}-\Sigma_{p,W}$. For convenience, we also include a
$K_{2}^{1/2}$ factor in $\|X\|_{W}$ to rescale the quantity to be of constant
order. For example, under our model (1.1) and some mild conditions, we show in
Section C of the Supplementary Material that $\|X\|_{W}=O(1)$ holds for a
general class of weight matrices.

In this paper, we define the cluster separation distance as
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$. First, we can view
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$ as a standardized distance, which
measures cluster separation per unit “variance” of $X$. Second,
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$ is invariant to the scale of $W$ and
also the scale of $X$ (e.g., transform $X$ to $tX$ for any $t\in\mathbb{R}$).
Finally, we note that our cluster separation metric depends on the choice of
the weight matrix $W$ since we use a weighted covariance distance to construct
the MCOD in (2.4). Thus, even if we consider the same data generating model
(1.1), the value of the cluster separation metric may differ depending on the
choice of $W$. This has important implications for clustering consistency (see
Remark 3.2 and Section D in the Supplementary Material) and the minimax lower
bound (see Section 3.2).

If the weight $W$ is known, we can directly apply our Algorithm 1 with the
input $\widehat{\Sigma}_{p,W}$. However, if $W$ depends on unknown parameters
(e.g., the weight $W_{O}=B(B^{T}B)^{-2}B^{T}/K_{2}$), we need to estimate $W$
and use a data dependent weight in Algorithm 1. To study clustering
consistency of the algorithm with a data dependent weight, we assume that
there exists an initial estimator $\widehat{W}$ of a deterministic weight
matrix $W$. To simplify the theoretical analysis, we focus on analyzing
Algorithm 1 with sample splitting. Specifically, we randomly divide the data
into two folds, $\\{X^{(i)}:i\in D_{1}\\}$ and $\\{X^{(i)}:i\in D_{2}\\}$,
where $D_{1}\cap D_{2}=\emptyset$ and $D_{1}\cup D_{2}=\\{1,...,n\\}$. The
estimator $\widehat{W}$ is constructed using the data in $D_{1}$, and then we
apply Algorithm 1 with the input
$\widehat{\Sigma}_{p,\widehat{W}}=\frac{1}{|D_{2}|}\sum_{i\in
D_{2}}X^{(i)}\widehat{W}X^{(i)T}$, where $\widehat{\Sigma}_{p,\widehat{W}}$ is
the weighted sample covariance matrix using the data in $D_{2}$. By using this
simple procedure, we can remove the dependence of $\widehat{W}$ and the data
in $D_{2}$.

Let us denote
$\Sigma_{p,\widehat{W}}=\mathbb{E}(X\widehat{W}X^{T}|\widehat{W})$, where the
expectation is taken with respect to $X$ and is independent of $\widehat{W}$.
The following main theorem in this section shows the clustering consistency of
our algorithm in the non-asymptotic regime.

######  Theorem 3.1.

Under model (1.1), assume that $\textrm{vec}(X)$ is multivariate Gaussian,
$\log p=o(n)$ and the following two conditions hold:

  * (A1)

Cluster separation condition:
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}>c_{0}\eta$, where $\eta\geq
c_{1}\sqrt{\frac{\log p}{nK_{2}}}$ for a universal constant $c_{1}>0$ and an
arbitrary constant $c_{0}\geq 4$.

  * (A2)

Stability condition:
$\Big{\\{}\mathrm{MCOD}(\Sigma_{p,\widehat{W}})/\|X\|_{\widehat{W}}\Big{\\}}\Big{/}\Big{\\{}\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}\Big{\\}}>4/c_{0}$,
where $c_{0}$ is defined in (A1).

Then using our Algorithm 1 with $\widehat{\Sigma}_{p,\widehat{W}}$ and the
threshold $\alpha=2\eta\cdot\|X\|_{\widehat{W}}$, we obtain perfect row
cluster recovery (i.e., $\widehat{\mathcal{G}}^{(r)}=\mathcal{G}^{(r)}$) with
probability greater than $1-\frac{c_{2}}{p}$ for some constant $c_{2}>0$.

In the following, we start from the discussion on the conditions in Theorem
3.1. The Gaussian assumption is imposed in order to derive a sharp bound for
$(\widehat{\Sigma}_{p,\widehat{W}}-\Sigma_{p,\widehat{W}})_{jk}$ when applying
the Hanson-Wright inequality, as we can decorrelate two dependent Gaussian
variables to make them independent. Recently, Adamczak (2015) derived a
variant of the Hanson-Wright inequality for dependent data with the so-called
convex concentration property. Using this new inequality, we can relax the
Gaussian assumption to more general distributions (e.g., sub-Gaussian) with
the convex concentration property. The conclusion in Theorem 3.1 remains the
same. However, to keep our presentation focused, we impose the Gaussian
assumption in this theorem.

The condition $\log p=o(n)$ is standard for high-dimensional data. We further
assume two major conditions, (A1) and (A2). Recall from the previous
discussion that the cluster separation is measured by
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$. Condition (A1) implies that in the
ideal case (i.e., the weight matrix $W$ is known), the clusters must be
separated by a factor of $\sqrt{\frac{\log p}{nK_{2}}}$. For this reason, we
call (A1) the cluster separation condition. On top of (A1), we also need
condition (A2), because in some cases the target weight matrix $W$ (e.g., the
optimal weight $W_{O}$) needs to be estimated. Condition (A2) quantifies the
stability of the cluster separation metric
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$ with respect to the perturbation of
$W$. Essentially, condition (A2) guarantees that
$\mathrm{MCOD}(\Sigma_{p,\widehat{W}})/\|X\|_{\widehat{W}}$ with the estimate
$\widehat{W}$ is still beyond the order $\sqrt{\frac{\log p}{nK_{2}}}$. When
$W$ does not need to be estimated, we can simply use $\widehat{\Sigma}_{p,W}$
in Algorithm 1, and in this case (A2) holds trivially. It is important to note
that the arbitrary constant $c_{0}$ in condition (A2) is $\geq$ 4, which
allows the separation based on $\widehat{W}$ to be smaller than that on $W$.
Indeed, the interplay between the two conditions (A1) and (A2) is
characterized by $c_{0}$. A larger $c_{0}$ can relax the stability condition
(A2), whereas the cluster separation condition (A1) becomes more stringent.
Lastly, we note that (A1) and (A2) are both high-level technical conditions,
which will be further explored in Section 4 under additional modeling
assumptions.

One important implication of Theorem 3.1 is that the minimum cluster
separation for clustering consistency is of order $\sqrt{\frac{\log
p}{nK_{2}}}$, which decreases as the number of column clusters $K_{2}$ grows.
In other words, clustering the rows of $X$ becomes easier if the columns of
$X$ have more column clusters. Indeed, this phenomenon is reasonable as the
data from two different column clusters show weaker dependence and therefore
improve the convergence rate of
$(\widehat{\Sigma}_{p,\widehat{W}}-\Sigma_{p,\widehat{W}})_{jk}$ in the
Hanson-Wright inequality. This result clearly demonstrates the benefit of
clustering the matrix $X$ over vector clustering. Finally, we note that in the
special case $q=K_{2}=1$, the order of our cluster separation metric matches
the existing result for vector clustering in Bunea et al. (2020).

Remark 3.2. The results in Theorem 3.1 are generally applicable to our
clustering algorithm with any positive semi-definite matrix $W$, provided (A1)
and (A2) hold. Recall that in Section 2.2, we consider two specific weights,
$W_{I}=I_{q}/q$ and $W_{O}=B(B^{T}B)^{-2}B^{T}/K_{2}$, where the latter was
called the optimal weight. We can apply our algorithm with either
$\widehat{\Sigma}_{p,W_{I}}$ or $\widehat{\Sigma}_{p,\widehat{W}_{O}}$, where
$\widehat{W}_{O}$ is an estimate of $W_{O}$ defined in Section 2.2. Since the
rate of the cluster separation for consistency is always $\sqrt{\frac{\log
p}{nK_{2}}}$, which does not depend on $W$, one might be tempted to conclude
that there is no benefit of using the optimal weight $W_{O}$ over $W_{I}$ for
clustering consistency. However, this conclusion is imprecise, as the value of
the cluster separation metric $\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$ depends
on $W$ and may differ substantially. When $q\gg K_{2}$, our algorithm with the
optimal weight attains clustering consistency in the presence of a larger
noise level $\sigma^{2}$ compared to when using $W_{I}$. A more detailed
illustration of this point can be found in Section D in the Supplementary
Material.

###  3.2 Minimax Optimality

In this section, we establish the minimax lower bound for the matrix
clustering problem. In this paper, we will focus on the optimal weight $W_{O}$
and from it construct an appropriate parameter space. Assume that
$\textrm{vec}(X)\sim N(0,\Sigma)$, where $\Sigma\in\mathbb{R}^{pq\times pq}$.
We define the relevant parameter space as

$M_{O}(p,q,K_{1},K_{2},\eta)\\\ =\big{\\{}\Sigma\in\mathbb{R}^{pq\times
pq}\big{|}\leavevmode\nobreak\ X\leavevmode\nobreak\ \textrm{satisfies
model}\leavevmode\nobreak\ (\ref{eq_model}),\leavevmode\nobreak\
\leavevmode\nobreak\
\mathrm{MCOD}(\Sigma_{p,W_{O}})/\|X\|_{W_{O}}\geq\eta\big{\\}},$

where the cluster separation metric is defined based on $W_{O}$. The following
theorem provides the lower bound for clustering over the parameter space
$M_{O}:=M_{O}(p,q,K_{1},K_{2},\eta)$.

######  Theorem 3.2.

For $K_{1}\geq 3$, there exists a positive constant $c$ such that, for any
$\eta$

$\displaystyle 0\leavevmode\nobreak\ \leq\leavevmode\nobreak\
\eta\leavevmode\nobreak\ <\leavevmode\nobreak\ c\sqrt{\frac{\log
p}{nK_{2}}}\leavevmode\nobreak\ ,$

we have

$\displaystyle\underset{\widehat{\mathcal{G}}}{\inf}\leavevmode\nobreak\
\underset{\Sigma\in M_{O}}{\sup}\leavevmode\nobreak\
\mathbb{P}_{\Sigma}(\widehat{\mathcal{G}}\neq\mathcal{G})\leavevmode\nobreak\
\geq\leavevmode\nobreak\ \frac{1}{7},$

where the infimum is taken over all possible estimators of $\mathcal{G}$.

Theorem 3.2 shows that it is impossible to attain clustering consistency
uniformly over the parameter space $M_{O}$ when the minimum cluster separation
value $\eta$ is below the threshold $\big{(}\frac{\log
p}{nK_{2}}\big{)}^{1/2}$, which matches the rate of $\eta$ in condition (A1)
in Theorem 3.1. Thus, provided the stability condition (A2) in Theorem 3.1
holds, our Algorithm 1 using $\widehat{\Sigma}_{p,\widehat{W}_{O}}$ and
$\alpha=2\eta\cdot\|X\|_{\widehat{W}_{O}}$ is minimax optimal for clustering.

We also have a minimax lower bound result in which the cluster separation
metric is defined with an arbitrary column membership matrix $\bar{B}$. The
setup of the appropriate parameter space, the presentation of the theorem, its
derivation and discussion can be found in Section E in the Supplementary
Material.

##  4 Applications to Matrix Normal Models

As seen from the previous section, our Algorithm 1 using the weighted
covariance matrix $\widehat{\Sigma}_{p,\widehat{W}_{O}}$ is minimax optimal,
provided conditions (A1) and (A2) in Theorem 3.1 hold. In this section, we
will verify conditions (A1) and (A2) with $\widehat{W}=\widehat{W}_{O}$. To
make the analysis of the cluster separation metric
$\mathrm{MCOD}(\Sigma_{p,W})/\|X\|_{W}$ tractable, we will make some
additional modeling assumptions.

On top of our model (1.1), we further assume that the latent variable $Z$
follows the matrix normal distribution. That is $Z\sim\mathrm{MN}(0,U,V)$,
where $U\in\mathbb{R}^{K_{1}\times K_{1}}$ and $V\in\mathbb{R}^{K_{2}\times
K_{2}}$ are symmetric positive definite matrices. Note that this is equivalent
to saying that $\mathrm{vec}(Z)\sim\mathrm{MVN}(\mathrm{vec}(0),V\otimes U)$,
which gives us
$\mathbb{E}(Z_{jk}Z_{j^{\prime}k^{\prime}})=U_{jj^{\prime}}V_{kk^{\prime}}$
and $\mathbb{E}(Z_{jk})=0$ for any $1\leq j,j^{\prime}\leq K_{1}$ and $1\leq
k,k^{\prime}\leq K_{2}$.

Recall that
$\widehat{W}_{O}=\widehat{B}(\widehat{B}^{T}\widehat{B})^{-2}\widehat{B}^{T}/s$,
where $\widehat{B}$ is an initial estimator of the column membership matrix
$B$ and $s$ is the estimated number of clusters. To facilitate the analysis,
we use $\mathcal{G}^{(c)}=\\{G_{1}^{(c)},...,G_{K_{2}}^{(c)}\\}$ and
$\widehat{\mathcal{G}}^{(c)}=\\{\widehat{G}_{1}^{(c)},...,\widehat{G}_{s}^{(c)}\\}$
to denote the true column cluster structure and the estimated column
structure, respectively. We define a $K_{2}\times s$ matrix $G$ that carries
information about the clustering accuracy of the initial estimator
$\widehat{B}$:

$G=B^{T}\widehat{B}(\widehat{B}^{T}\widehat{B})^{-1}=\begin{bmatrix}\frac{|G_{1}^{(c)}\cap\widehat{G}_{1}^{(c)}|}{|\widehat{G}_{1}^{(c)}|}&\frac{|G_{1}^{(c)}\cap\widehat{G}_{2}^{(c)}|}{|\widehat{G}_{2}^{(c)}|}&...&\frac{|G_{1}^{(c)}\cap\widehat{G}_{s}^{(c)}|}{|\widehat{G}_{s}^{(c)}|}\\\
:&:&&:\\\
\frac{|G_{K_{2}}^{(c)}\cap\widehat{G}_{1}^{(c)}|}{|\widehat{G}_{1}^{(c)}|}&\frac{|G_{K_{2}}^{(c)}\cap\widehat{G}_{2}^{(c)}|}{|\widehat{G}_{2}^{(c)}|}&...&\frac{|G_{K_{2}}^{(c)}\cap\widehat{G}_{s}^{(c)}|}{|\widehat{G}_{s}^{(c)}|}\\\
\end{bmatrix}.$

Note that the columns of $G$ sum to $1$, and in the ideal case of
$\mathcal{G}^{(c)}=\widehat{\mathcal{G}}^{(c)}\text{ (i.e., }\widehat{B}=B)$,
we have that $G=I_{K_{2}}$. Denote by

$C_{K}=\max_{1\leq a\leq
p}\frac{1}{K_{2}}\sum_{t=1}^{K_{2}}\frac{(\sum_{j\in[t]}\sigma^{2}_{aj})^{2}}{|[t]|^{4}}$
(4.1)

the maximum weighted average of the error variances over $K_{2}$ column
clusters, where $[t]$ denotes the $t$-th column cluster. Similarly, we define

$C_{s}=\max_{1\leq a\leq
p}\frac{1}{s}\sum_{t=1}^{s}\frac{(\sum_{j\in[\widehat{t}]}\sigma^{2}_{aj})^{2}}{|[\widehat{t}]|^{4}}$
(4.2)

where $[\widehat{t}]$ denotes the estimated $t$-th column cluster. The
following proposition shows the conditions under which (A1) and (A2) in
Theorem 3.1 hold in the matrix normal model. The proof can be found in Section
H.3.1 of the Supplementary Material.

######  Proposition 4.1.

Under model (1.1) and the above matrix normal model assumptions, if we assume:

  1. (P0)

$C_{\min}\leavevmode\nobreak\ \leq\leavevmode\nobreak\
\lambda_{\min}(V)\leavevmode\nobreak\ \leq\leavevmode\nobreak\
\lambda_{\max}(V)\leavevmode\nobreak\ \leq\leavevmode\nobreak\ C_{\max}$ for
two positive constants $C_{\min}$ and $C_{\max}$.

  2. (P1)

(Cluster Separation) For any $1\leq j\neq k\leq K_{1}$, there exists $1\leq
l\leq K_{1}$ such that

$\displaystyle\big{|}U_{jl}-U_{kl}\big{|}\leavevmode\nobreak\ $
$\displaystyle>\leavevmode\nobreak\ c_{0}\cdot
c_{1}\sqrt{\frac{\log(p)}{nK_{2}}}\cdot\frac{K_{2}}{\text{tr}(V)}\cdot\Bigg{\\{}||\text{diag}(U)||_{\max}\cdot
C_{\max}+C_{K}^{1/2}\Bigg{\\}}$

holds for a universal constant $c_{1}$ and an arbitrary constant $c_{0}\geq
4$, where
$||\text{diag}(U)||_{\max}=\max_{1\leq j\leq K_{1}}U_{jj}$.

  3. (P2)

(Stability) Either (i) or (ii) holds:

$\text{(i)}\hskip
28.45274pt\begin{cases}\sqrt{\min(s,K_{2})}\cdot|\lambda_{\max}(GG^{T})|\cdot||\text{diag}(U)||_{\max}\cdot
C_{\max}\leavevmode\nobreak\ \leq\leavevmode\nobreak\
\sqrt{s}\cdot\sqrt{C_{s}}\\\
\frac{1}{\lambda_{\min}(GG^{T})}\leavevmode\nobreak\ \leq\leavevmode\nobreak\
\frac{c_{0}}{8}\cdot\sqrt{\frac{C_{K}}{C_{s}}}\cdot\sqrt{\frac{K_{2}}{s}}\end{cases}$
$\text{(ii)}\hskip
28.45274pt\begin{cases}\sqrt{\min(s,K_{2})}\cdot|\lambda_{\max}(GG^{T})|\cdot||\text{diag}(U)||_{\max}\cdot
C_{\max}\leavevmode\nobreak\ >\leavevmode\nobreak\
\sqrt{s}\cdot\sqrt{C_{s}}\\\
\frac{\lambda_{\max}(GG^{T})}{\lambda_{\min}(GG^{T})}\leavevmode\nobreak\
\leq\leavevmode\nobreak\
\frac{c_{0}}{8}\cdot\frac{C_{\min}}{C_{\max}}\cdot\frac{||\text{diag}(U)||_{\min}}{||\text{diag}(U)||_{\max}}\cdot\sqrt{\frac{K_{2}}{\min(s,K_{2})}}\end{cases}$

where $c_{0}$ is the arbitrary constant from (P1), $C_{K}$ and $C_{s}$ are
defined in (4.1) and (4.2), and $||\text{diag}(U)||_{\min}=\min_{1\leq j\leq
K_{1}}U_{jj}$.

Then the conditions (A1) and (A2) in Theorem 3.1 hold.

The assumption (P0) is standard and commonly used in the high-dimensional
statistics literature (Basu and Michailidis, 2015; Cai et al., 2010; Bai and
Silverstein, 2010). For (P1), it says that any two rows of the row covariance
matrix $U$ cannot be nearly identical because otherwise the corresponding two
clusters would not be identifiable. When the error variances are not too large
(e.g., $\sum_{j\in[t]}\sigma^{2}_{aj}\lesssim|[t]|^{2}$), then $C_{K}$ behaves
like a constant. If we further assume $||\text{diag}(U)||_{\max}=O(1)$, then
the cluster separation reduces to
$\big{|}U_{jl}-U_{kl}\big{|}\gtrsim\sqrt{\frac{\log(p)}{nK_{2}}}$, which
decreases with $n$ and $K_{2}$ as we have discussed before. Finally, (P2)
requires the estimated cluster structure to be close enough to the true
cluster structure. When the initial estimator $\widehat{B}$ is reasonably
accurate, we would expect $C_{K}$ and $C_{s}$ to be of the same order. When
$K_{2}\asymp s$ and
$||\text{diag}(U)||_{\max}\asymp||\text{diag}(U)||_{\min}$, (P2) reduces to
the condition that the smallest eigenvalue of $GG^{T}$ is bounded from below
by a constant and its largest eigenvalue is bounded from above by a constant.
Recall that when $\widehat{B}$ is more accurate, $G$ would be closer to an
identity matrix and therefore (P2) would be more likely to hold. Thus, (P2)
essentially gives a sharp characterization of the “contraction region” in
which the algorithm with the imperfect initial estimator $\widehat{B}$ still
leads to clustering consistency. In fact, the following proposition gives a
concrete example where $\widehat{B}$ satisfies the condition (P2).

######  Proposition 4.2.

Under the matrix normal model, assume that the noise variances are homogeneous
and the cluster sizes are balanced, (i.e., $\forall i,j$,
$\mathbb{E}(\Gamma_{ij})=\sigma^{2}$, $\frac{M_{q}}{m_{q}}\leq C$ for some
constant $C\geq 1$ where $M_{q}$ and $m_{q}$ are the largest and smallest
column cluster sizes, respectively). If $c_{0}\geq
8C\big{/}\big{(}\frac{C_{\min}}{C_{\max}}\cdot\frac{||\text{diag}(U)||_{\min}}{||\text{diag}(U)||_{\max}}\big{)}$
and $\sqrt{M_{q}}\geq\frac{\sigma^{2}}{||\text{diag}(U)||_{\max}\cdot
C_{\max}}$, then the stability condition (P2) in Proposition 4.1 is satisfied
with $\widehat{B}=I_{q}$.

This initial estimate $\widehat{B}=I_{q}$, which in turn gives us
$\leavevmode\nobreak\ \widehat{W}_{O}\leavevmode\nobreak\
=\leavevmode\nobreak\
\frac{1}{q}\bar{B}(\bar{B}^{T}\bar{B})^{-2}\bar{B}^{T}\leavevmode\nobreak\
=\leavevmode\nobreak\ \frac{1}{q}I_{q}$ is actually equivalent to the $W_{I}$
that was considered in Section 2.2. This “naive” weight that assigns each
column to its own column cluster is practical, since it does not require any
prior knowledge on the column cluster structure. This is significant as it
guarantees a feasible minimax optimal implementation of our Algorithm 1 for
matrix normal models. The proof can be found in Section H.3.2 of the
Supplementary Material.

##  5 Some Practical Considerations

In this section, we discuss the practical implementation of our algorithm
using the optimal weight. First, we look at the iterative variant of our
algorithm and afterwards we discuss miscellaneous aspects of the algorithm
such as sample splitting, standardization and selecting the tuning parameter
$\alpha$.

###  5.1 The Iterative One-Step and Two-Step Methods

We first introduce the iterative one-step hierarchical algorithm with the
optimal weight in Algorithm 2. Skipping over the details on sample splitting
which were already discussed in Section 3.1 and will be further discussed in
Section 5.2, we first apply our hierarchical Algorithm 1 to cluster the rows
using $\bar{B}=I_{q}$ as the initial estimate for $B$ and
$\widehat{W}_{O,(r)}=\frac{1}{q}I_{q}$ as an estimate of the optimal weight
matrix for the rows. This is the same as the initial estimate proposed in
Proposition 4.2. Then we use the obtained $\widehat{A}$ to construct an
estimate of the optimal weight matrix for the columns
$\widehat{W}_{O,(c)}=\widehat{A}(\widehat{A}^{T}\widehat{A})^{-2}\widehat{A}^{T}/t$
(where $t$ is the number of clusters from $\widehat{A}$), and use it in
Algorithm 1 to cluster the rows to get $\widehat{B}$. We provide the
theoretical guarantees on consistency and minimax optimality for the iterative
one-step Algorithm 2 in Section F.1.1 and Section F.1.2, respectively, of the
Supplementary Material.

Intuitively, we can further repeat steps 2-3 in Algorithm 2 to iteratively
cluster the rows and columns of $X$, leading to a multi-step algorithm. The
two-step algorithm essentially repeats step 2 in Algorithm 2 one more time
after step 3. The details can be found in Section F.2 of the Supplementary
Material.

  1. 1.

Split the data into two folds: $D_{1}$, $D_{2}$. Set an initial value
$\bar{B}$ (e.g., $\bar{B}=I_{q}$).

  2. 2.

On data $D_{1}$, apply Algorithm 1 with
$\widehat{\Sigma}_{p,\widehat{W}_{O,(r)}}=\frac{1}{|D_{2}|}\sum_{i\in
D_{2}}X^{(i)}\widehat{W}_{O,(r)}X^{(i)T}$ to

cluster the rows of $X$, where
$\widehat{W}_{O,(r)}=\bar{B}(\bar{B}^{T}\bar{B})^{-2}\bar{B}^{T}/s$ and $s$
denotes the

estimated number of clusters from $\bar{B}$. Obtain the resulting row cluster
$\widehat{\mathcal{G}}^{(r)}$ or equivalently the membership matrix
$\widehat{A}$.

* 3.

On data $D_{2}$, apply Algorithm 1 with
$\widehat{\Sigma}_{q,\widehat{W}_{O,(c)}}=\frac{1}{|D_{3}|}\sum_{i\in
D_{3}}X^{(i)T}\widehat{W}_{O,(c)}X^{(i)}$ to

cluster the columns of $X$, where
$\widehat{W}_{O,(c)}=\widehat{A}(\widehat{A}^{T}\widehat{A})^{-2}\widehat{A}^{T}/t$
and $t$ denotes the

estimated number of clusters from $\widehat{A}$. Obtain the resulting cluster
$\widehat{\mathcal{G}}^{(c)}$ or

equivalently the membership matrix $\widehat{B}$.

Algorithm 2 Iterative One-Step Hierarchical Algorithm with Optimal Weight

###  5.2 Sample Splitting, Standardization and Selecting $\alpha$

We note that in Section 3.1, to facilitate the theoretical analysis of our
algorithm, we split the data into two folds. Such a sample splitting procedure
is feasible in practice when $n$ is relatively large. However, in practice,
when we apply our algorithm to the data with very small $n$ (e.g., $n=40$),
sample splitting tends to yield unstable clustering results. Considering this,
in practice, we chose to implement both steps 2 and 3 in Algorithm 2 on the
entire dataset. A detailed comparison of the performance of our method with
and without sample splitting is presented in a simulation study in Section F.3
of the Supplementary Material.

Another point to mention is that in practice, the feature matrix $X$ may have
elements with differing variances. While the theoretical guarantees of our
algorithm in Section 3.1 remain valid, the concentration bound for
$(\widehat{\Sigma}_{p,W}-\Sigma_{p,W})_{jk}$ via the Hanson-Wright inequality
will be dominated by the variables with a larger variance. To tighten this
upper bound, we recommend applying our algorithm to standardized data - $X$
whose elements all have mean 0 and variance 1. Empirically, we observe that
clustering accuracy can be significantly improved when the algorithm is
applied to the standardized data.

Finally, we note that, in steps 2 and 3 in Algorithm 2, we need to choose the
threshold value $\alpha$ in our hierarchical algorithm. In Section F.4 of the
Supplementary Material, we present Algorithm 6, a data-driven cross-validation
scheme to choose the optimal threshold value $\alpha$.

##  6 Further Extensions

In this section we discuss further extensions of our model (1.1) and COD based
methods of Algorithm 1, 2 and 3 to account for a variety of more general
settings.

###  6.1 The Dependent Noise Model

In model (1.1), we assume that the entries in the noise matrix $\Gamma$ are
independent. In practice, the model can be more flexible by allowing the
entries of $\Gamma$ to have a dependence structure. In this section, we
consider the more general model with correlated noise variables and present
the conditions needed for clustering consistency in Theorem 6.1. Define

$\displaystyle\gamma(\Sigma_{p,W})\leavevmode\nobreak\ \leavevmode\nobreak\
:=\leavevmode\nobreak\ \leavevmode\nobreak\ \underset{1\leq a,b\leq
p}{\max}\leavevmode\nobreak\ \leavevmode\nobreak\ \underset{c\neq
a,b}{\max}\leavevmode\nobreak\ \leavevmode\nobreak\
\Big{|}\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{ac}-\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{bc}\Big{|}.$ (6.1)

This quantity $\gamma(\Sigma_{p,W})$ plays an important role in determining
how difficult the clustering problem is when using the COD method with
correlated noise variables. More specifically, it contains information on the
maximum difference of two different elements in the same row or column in the
$p\times p$ weighted covariance matrix. Note that in the independent noise
setting, $\gamma(\Sigma_{p,W})=0$ holds since all off-diagonal terms of
$\mathbb{E}(\Gamma W\Gamma^{T})$ would be 0. The analog to the clustering
consistency result in Theorem 3.1 under this more general noise setting using
$\gamma(\Sigma_{p,W})$ is outlined in the following theorem.

######  Theorem 6.1.

(Consistency with Dependent Noise Variables)
Under the model $X=AZB^{T}+\Gamma$, assume that $\textrm{vec}(X)$ is
multivariate Gaussian, $\log p=o(n)$ and the following two conditions hold:

  * (A1)

Cluster separation condition:

$\text{MCOD}(\Sigma_{p,W})/\|X\|_{W}\leavevmode\nobreak\ \leavevmode\nobreak\
>\leavevmode\nobreak\ \leavevmode\nobreak\ c_{0}\eta\leavevmode\nobreak\
+\leavevmode\nobreak\ \frac{\gamma(\Sigma_{p,W})}{||X||_{W}},$ (6.2)

where $c_{0}\geq 4$ is an arbitrary constant, $\leavevmode\nobreak\ \eta\geq
c_{1}\sqrt{\frac{\log p}{nK_{2}}}\leavevmode\nobreak\ $ for a universal
constant $c_{1}>0$ and $\gamma(\Sigma_{p,W})$ is defined in (6.1).

  * (A2)

Stability condition:

$\bigg{[}\Big{\\{}\text{MCOD}(\Sigma_{p,\widehat{W}})-\gamma(\Sigma_{p,\widehat{W}})\Big{\\}}\Big{/}\|X\|_{\widehat{W}}\bigg{]}\Bigg{/}\bigg{[}\Big{\\{}\text{MCOD}(\Sigma_{p,W})-\gamma(\Sigma_{p,W})\Big{\\}}\Big{/}\|X\|_{W}\bigg{]}>\frac{4}{c_{0}}$

where $c_{0}$ is defined in (A1).

Then using our Algorithm 1 with $\widehat{\Sigma}_{p,\widehat{W}}$ and the
threshold $\alpha=2\eta\cdot\|X\|_{\widehat{W}}\leavevmode\nobreak\
+\leavevmode\nobreak\ \gamma(\Sigma_{p,\widehat{W}})$, we obtain perfect row
cluster recovery (i.e., $\widehat{\mathcal{G}}^{(r)}=\mathcal{G}^{(r)}$) with
probability greater than $1-\frac{c_{2}}{p}$ for some constant $c_{2}>0$.

Compared to Theorem 3.1, the cluster separation is stricter as the
$\frac{\gamma(\Sigma_{p,W})}{||X||_{W}}$ factor on the right hand side of
(6.2) shows that MCOD has to be much larger to account for the correlated
noise structure. Further discussions of this result can be found in Section
G.1 of the Supplementary Material.

###  6.2 Nested Clustering to Incorporate Mean Information

Our proposed model (1.1), Algorithm 1, 2 and 3 on their own can have many uses
in practice where the dependence structure of the data is of main importance.
However, there may be instances where mean information is also informative. In
order to take full advantage of both components, we generalize our model and
propose a method to accommodate both the mean and covariance information in
clustering. More specifically, we propose a non-centered latent variable model
$X=M+AZB^{T}+\Gamma$, where $M=\mathbb{E}(X)$ induces the first layer row and
column cluster structures, and the membership matrices $A$ and $B$ encode the
second layer row and column cluster structures (from the covariance of $X$),
which is assumed to be nested inside the mean clusters. We further propose a
nested clustering algorithm in which a mean-based clustering method can be
implemented first, and then on each cluster, our covariance-based method can
be applied to capture the finer, more intricate partitions within each broad
cluster. We also include a simulation study comparing the empirical
performance of this nested method and our original Algorithm 3. The full
presentation can be found in Section G.2 of the Supplementary Material.

###  6.3 Higher Order Tensor Models

Recall that our original problem is to cluster $n$ i.i.d. $p\times q$ matrices
over the $p$ and $q$ directions. We now discuss the extension our model to
higher order tensor settings. Consider a three-way tensor
$X\in\mathbb{R}^{J\times P\times Q}$, which satisfies
$X=Z\times_{1}A\times_{2}B\times_{3}C\leavevmode\nobreak\
+\leavevmode\nobreak\ \Gamma$, where $Z\in\mathbb{R}^{K_{1}\times K_{2}\times
K_{3}}$ is a three-way latent tensor with $\mathbb{E}(Z)=0$,
$A\in\mathbb{R}^{J\times K_{1}}$, $B\in\mathbb{R}^{P\times K_{2}}$ and
$C\in\mathbb{R}^{Q\times K_{3}}$ are the unknown membership matrices for the
rows, columns and tubes of $X$, respectively, and
$\Gamma\in\mathbb{R}^{J\times P\times Q}$ represents the mean 0 random noise
tensor. Here, $Z\times_{1}A$ is the tensor $n$-mode product defined in Section
G.3 of the Supplementary Material. Given $n$ i.i.d copies of $X$, we
generalize our algorithm to recover the membership matrices $A,B$ and $C$. We
also conduct a simulation study to assess the empirical performance of our
algorithm. The full presentation can be found in Section G.3 of the
Supplementary Material.

##  7 Simulation Results

We consider the following data generating process. We fix $p=q=100$,
$K_{1}=K_{2}=10$ and a moderately unbalanced row and column cluster size
structure (both having cluster sizes of $3$, $6$, $6$, $8$, $10$, $10$, $12$,
$12$, $14$, $19$, respectively). This gives us the membership matrices $A$ and
$B$. The latent variable $Z$ is generated from the matrix normal distribution
$\mathrm{MN}(0_{K_{1}\times K_{2}},U_{K_{1}\times K_{1}},V_{K_{2}\times
K_{2}})$, where $U_{jk}=(-0.4)^{|j-k|}$ and $V_{jk}=0.3^{|j-k|}$. We further
generate $\Gamma_{ij}\sim N(0,\sigma_{ij}^{2})$, where we consider the
following three settings for the noise variance: (1) homogeneous noise
variances $\sigma_{ij}^{2}=15$; (2) heterogeneous noise variances proportional
to the corresponding row and column cluster sizes:
$\sigma_{ij}^{2}=\frac{15pq\cdot
v_{ij}}{\sum_{i,j}v_{ij}},v_{ij}=\frac{m_{p}^{(i)}\cdot
m_{q}^{(j)}}{\sqrt{\frac{pq}{K_{1}K_{2}}}}$; and (3) heterogeneous noise
variances randomly generated from the Uniform distribution:
$\sigma_{ij}^{2}=\frac{15pq\cdot
u_{ij}^{h}}{\sum_{i,j}u_{ij}^{h}},u_{ij}\sim\mathrm{Unif}(0,1)$ where $h$
determines the level of heterogeneity. They will be referred to as the
“homogeneous”, “proportional” and “random” cases, respectively.

Recall that in Section 3.1, the noise variances play an important role on the
clustering consistency of our algorithm. We design these three settings to
test the performance of our algorithm under different patterns of noise
variances. We keep the mean noise variance to be 15 for all three cases. For
the third setting, we set $h=0.87$ because it gives us a similar level of
heterogeneity as the second setting. In both the second and third settings,
the standard deviation of the noise variances is around 7.95. Finally, we
generate $X$ from the model (1.1). We vary the sample size $n$ in the
simulations and the simulations are repeated 30 times. To measure clustering
accuracy, we consider the adjusted rand index (ARI) (Hubert and Arabie, 1985).
Note that an ARI value of 1 implies a perfect match between the true and
estimated cluster partitions. The formal definition of the ARI is shown in
Section K of the Supplementary Material.

We consider the following clustering methods: (1) DEEM; (2) TGMM; (3) TEMM;
(4) Algorithm 1 with the initial weight $W_{I}$ (NAIVE COD); (5) Algorithm 2
(1-STEP COD); and (6) Algorithm 3 (2-STEP COD in Section F.2 of the
Supplementary Material). We note that the data-driven cross-validation scheme
is used to choose the tuning parameter in our COD-based algorithms. The
competing methods we used in our simulations came from recent model based
tensor clustering papers, Mai et al. (2022) and Deng and Zhang (2022). Their
DEEM and TGMM methods utilize the tensor normal mixture model and their TEMM
method utilizes the tensor envelope mixture model. Together with some variants
of the EM algorithm, the authors exploited not only the mean information, but
also the covariance structure in each direction of the tensor. These methods
were implemented in our simulations using their designated functions within
the R package TensorClustering.

Since the goal of DEEM, TGMM and TEMM is to cluster over $n$ i.i.d copies of
tensors rather than the features, we need to adapt their methods for row and
column clustering in our simulations. In particular, we implement row
clustering by regarding the $p$ rows in our setting to be the direction of the
observations in their tensor setting, and regard our $q\times n$ matrix slices
to be the tensors in their tensor setting. Note that this violates their
tensor normal distributional assumptions as the $p$ “copies” of $q\times n$
matrix slices are no longer i.i.d. The column clustering is performed
similarly. The simulation results for the random noise variance setting are
shown in Figure 1. The results for the other two noise variance settings can
be found in Section I.1 of the Supplementary Material.

It must be noted that since the distributional assumptions for the DEEM, TGMM
and TEMM methods are violated, these methods fail to converge in some of the
runs, especially when $n$ is large. One possible reason is that when $n$ is
larger, the i.i.d. assumption is more severely violated and the computation
becomes indeterminate in the process. In order to fairly compare our method
with DEEM, TGMM and TEMM, we report not only the average ARI among the
successful runs of their methods, but also the percentage of the successful
runs out of the total number of runs. The latter is called the “successful
implementation rate”.

The simulation results show that in all three noise variance settings the
successful implementation rates of DEEM, TGMM and TEMM decrease drastically as
$n$ increases. In contrast, our COD methods never fail to converge in any of
the settings. When the competing methods converge, the ARI is decent,
especially when $n\leq 12$ where they outperform the COD methods in the
homogeneous and random noise variance settings. However, as $n$ grows, our COD
methods become more reliable and more accurate.

In most cases, our iterative algorithms 1-STEP COD and 2-STEP COD improve the
performance of NAIVE COD, which is consistent with our theoretical analysis.
In particular, 1-STEP COD and 2-STEP COD can achieve an ARI value close to 1
when $n=18$, whereas the NAIVE COD may require $n$ much larger than $30$ to
attain the same level of accuracy. The phenomenon holds for both row and
column clustering. We also find that when $n$ is moderate (e.g., $n\geq 18$),
1-STEP COD and 2-STEP COD have very similar performances. Due to the extra
computational cost of 2-STEP COD, we generally recommend 1-STEP COD for
practical use if $n$ is moderate or large. However if $n$ is small, 2-STEP COD
may outperform 1-STEP COD.

We conducted additional simulation studies to evaluate the performance of our
clustering algorithm with an unbalanced matrix structure of $p=100$ and
$q=15$. The detailed results can be found in Section I.2 of the Supplementary
Material.

Figure 1: The ARI and the successful implementation rates of our hierarchical
clustering Algorithm 1 with $W_{I}$ (NAIVE COD), the one-step Algorithm 2
(1-STEP COD), the two-step Algorithm 3 (2-STEP COD), and competing model-based
tensor clustering methods DEEM [Mai et al. (2022)], TGMM and TEMM [Deng and
Zhang (2022)] under the random noise variance setting. Here, $p=q=100$,
$K_{1}=K_{2}=10$, and we have moderately unbalanced cluster sizes of
$(3,6,6,8,10,10,12,12,14,19)$ for the rows and the columns each. Successful
implementation rates are especially emphasized as the competing methods (DEEM,
TGMM, TEMM) tended to fail during execution for some of the runs. The ARI
rates for these methods were averaged over the successful runs, and when there
were no successes, the ARI was reported as 0 with a dashed line.

##  8 Genomic Data Analysis

We apply our method to the atlas of gene expression in the mouse aging project
dataset (Zahn et al., 2007), which contains gene expression values of 8932
genes in 16 tissues for 40 mice. Similar to Ning and Liu (2013) and Yin and Li
(2012), we only focus on a subset of the data belonging to the mouse vascular
endothelial growth factor signaling pathway. In order to maximize our usage of
the tissue data, we drop 4 samples that have missing data for the tissues. So,
in the end, our data is a $37\times 12\times 36$ array that corresponds to
$p=37$ genes, $q=12$ tissues, and $n=36$ mice.

Figure 2 shows the gene clusters and tissue clusters obtained from our two-
step Algorithm 3 with a data-driven tuning parameter. Disregarding the
singletons - clusters with only one element - the estimated number of gene
clusters is 4 and the estimated number of tissue clusters is 3. In Zahn et al.
(2007), the authors group the tissues through hierarchical clustering, but
they group genes that are similarly age-regulated through empirical meta
analysis. Our method, on the other hand, can conveniently be applied to both
the genes and the tissues at the same time. Interestingly, the tissue
clustering result from our method agrees with the tissue clustering result in
Zahn et al. (2007) which classifies tissues into 3 groups - vascular, neural,
and steroid responsive. The “Cerebrum”, “Cerebellum”, and “Hippocampus” are
all neural tissues that are parts of the brain, and they are all clustered
into one tissue cluster by our algorithm. Also, “Adrenal” and “Thymus” tissues
are steroid responsive tissues, while the “Lung” and “Kidney” tissues are
vascular tissues. The respective pairs appear in our clustering result as
well. As for gene clustering, it is known in the biology community that the
O03Rik gene inhibits the Nfat5 gene. More specifically, O03Rik acts as a
negative regulator of the calcineurin/NFAT signaling pathway and it inhibits
NFAT nuclear translocation and transcriptional activity by suppressing the
calcium-dependent calcineurin phosphatase activity (OriGene, 2020). It is then
reasonable that the two genes are clustered together, because they would be
highly negatively correlated. Also, calcineurin (Ppp3) is a serine/threonine
protein phosphatase that is dependent on calcium and calcium modulated
proteins (Hogan et al., 2003). It activates nuclear factor of activated T cell
cytoplasmic (Nfatc), a transcription factor, by dephosphorylating it. This
Nfatc happens to be encoded in part by the Nfatc4 gene. On the other hand,
there are three isozymes of the catalytic subunit of calcineurin and this is
encoded in part by the Ppp3r1 gene (Liu et al., 2005). Thus, the Ppp3r1 and
Nfatc4 genes are both related to calcineurin - one having to do with how it is
generated, and one having to do with something it activates. This reasonably
explains the highly correlated results obtained from our clustering algorithm.
Another possible explanation is given in Heit et al. (2006) which postulates
that calcineurin/NFAT signaling is critical in $\beta$-cell growth in the
pancreas. All this external evidence suggests that our clustering results are
biologically meaningful. We also have cluster results when applying a
competing method DEEM (Deng and Zhang, 2022) on the same data. The cluster
results along with a comparative discussion are presented in Section J of the
Supplementary Material.

Figure 2: Gene clusters and tissue clusters obtained from our Algorithm 3.

##  9 Discussion

In this paper, we study the variable clustering problem for matrix valued data
$X$ under the latent variable model (1.1). A class of hierarchical clustering
algorithms based on the weighted covariance difference is proposed. The
theoretical and empirical performance of the algorithm heavily depends on the
weight matrix $W$. Theoretically, we show that under mild conditions, our
algorithm with a large class of weight matrices $W$ can attain clustering
consistency with high probability. To further characterize the effect of the
weight matrix $W$ on the theoretical performance of our algorithm, we
establish the minimax lower bound under the latent variable model (1.1), from
which we prove that our algorithm using the weight matrix $W_{O}$ is minimax
optimal for variable clustering. In particular, we apply the theory to the
more concrete matrix normal model and show that clustering consistency and
minimax optimality can be achieved in practice with an implementable initial
weight,
$\widehat{W}_{O,(r)}=\frac{1}{q}\bar{B}(\bar{B}^{T}\bar{B})^{-2}\bar{B}^{T}=\frac{1}{q}I_{q}$.
Empirically, we develop iterative one-step and two-step algorithms based on
the weight matrix $W_{O}$, which outperform competing methods.

While we introduced extensions of our latent variable model and algorithms in
Section 6, it would be interesting to further investigate each topic on its
own. In particular, the theoretical effect of using different weight matrices
when clustering higher order tensors with COD based methods would be an
interesting direction. Another interesting extension would be to generalize
our method to the overlapping clustering setting (Bing et al., 2020), where
the rows and columns may simultaneously belong to multiple row and column
clusters. For now, however, we leave them for future study.

## Acknowledgment

Ning is supported in part by National Science Foundation (NSF) CAREER award
DMS-1941945, NSF award DMS-1854637 and NIH 1RF1AG077820-01A1.

## Supplementary Material

The supplement consists of additional technical details, derivations and
discussions. Relevant proofs and additional numerical results are also
included.

## References

  * Abbe (2017) Abbe, E. (2017).  Community detection and stochastic block models: recent developments.  The Journal of Machine Learning Research 18 6446–6531.
  * Adamczak (2015) Adamczak, R. (2015).  A note on the Hanson-Wright inequality for random vectors with dependencies.  Electronic Communications in Probability 20 1 – 13\.
  * Airoldi et al. (2008) Airoldi, E. M., Blei, D., Fienberg, S. and Xing, E. (2008).  Mixed membership stochastic blockmodels.  Advances in neural information processing systems 21.
  * Bai and Silverstein (2010) Bai, Z. and Silverstein, J. W. (2010).  Spectral analysis of large dimensional random matrices, vol. 20.  Springer.
  * Basu and Michailidis (2015) Basu, S. and Michailidis, G. (2015).  Regularized estimation in sparse high-dimensional time series models.  The Annals of Statistics 43 1535–1567.
  * Bing et al. (2020) Bing, X., Bunea, F., Ning, Y. and Wegkamp, M. (2020).  Adaptive estimation in structured factor models with applications to overlapping clustering.  The Annals of Statistics 48 2055–2081.
  * Bunea et al. (2020) Bunea, F., Giraud, C., Luo, X., Royer, M. and Verzelen, N. (2020).  Model assisted variable clustering: minimax-optimal recovery and algorithms.  Annals of statistics 48 111\.
  * Cai et al. (2010) Cai, T. T., Zhang, C.-H. and Zhou, H. H. (2010).  Optimal rates of convergence for covariance matrix estimation.  The Annals of Statistics 38 2118–2144.
  * Chi et al. (2017) Chi, E. C., Allen, G. I. and Baraniuk, R. G. (2017).  Convex biclustering.  Biometrics 73 10–19.
  * Chi et al. (2020) Chi, E. C., Gaines, B. R., Sun, W. W., Zhou, H. and Yang, J. (2020).  Provable convex co-clustering of tensors.  The Journal of Machine Learning Research 21 8792–8849.
  * Deng and Zhang (2022) Deng, K. and Zhang, X. (2022).  Tensor envelope mixture model for simultaneous clustering and multiway dimension reduction.  Biometrics 78 1067–1079.
  * Eisenach et al. (2020) Eisenach, C., Bunea, F., Ning, Y. and Dinicu, C. (2020).  High-dimensional inference for cluster-based graphical models.  Journal of Machine Learning Research 21 1–55.
  * Everitt et al. (2011) Everitt, B. S., Landau, S., Leese, M. and Stahl, D. (2011).  Cluster analysis 5th ed.  John Wiley.
  * Fraiman and Li (2022) Fraiman, N. and Li, Z. (2022).  Biclustering with alternating k-means.
URL https://doi.org/10.48550/arXiv.2009.04550

  * Hartigan (1972) Hartigan, J. A. (1972).  Direct clustering of a data matrix.  Journal of the american statistical association 67 123–129.
  * Hastie et al. (2009) Hastie, Tibshirani and Friedman (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition.  Springer Series in Statistics, Springer.
  * Heit et al. (2006) Heit, J. J., Apelqvist, Å. A., Gu, X., Winslow, M. M., Neilson, J. R., Crabtree, G. R. and Kim, S. K. (2006).  Calcineurin/nfat signalling regulates pancreatic $\beta$-cell growth and function.  Nature 443 345–349.
  * Hogan et al. (2003) Hogan, P. G., Chen, L., Nardone, J. and Rao, A. (2003).  Transcriptional regulation by calcium, calcineurin, and nfat.  Genes & development 17 2205–2232.
  * Holland et al. (1983) Holland, P. W., Laskey, K. B. and Leinhardt, S. (1983).  Stochastic blockmodels: First steps.  Social networks 5 109–137.
  * Hubert and Arabie (1985) Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of classification 2 193–218.
  * Kogan (2007) Kogan, J. (2007).  Introduction to clustering large and high-dimensional data.  Cambridge University Press.
  * Lee et al. (2010) Lee, M., Shen, H., Huang, J. Z. and Marron, J. (2010).  Biclustering via sparse singular value decomposition.  Biometrics 66 1087–1095.
  * Li et al. (2010) Li, X., Ye, Y., Li, M. J. and Ng, M. K. (2010).  On cluster tree for nested and multi-density data clustering.  Pattern Recognition 43 3130–3143.
  * Liu et al. (2005) Liu, L., Zhang, J., Yuan, J., Dang, Y., Yang, C., Chen, X., Xu, J. and Yu, L. (2005).  Characterization of a human regulatory subunit of protein phosphatase 3 gene (ppp3rl) expressed specifically in testis.  Molecular biology reports 32 41–45.
  * Madeira and Oliveira (2004) Madeira, S. C. and Oliveira, A. L. (2004).  Biclustering algorithms for biological data analysis: a survey.  IEEE/ACM transactions on computational biology and bioinformatics 1 24–45.
  * Mai et al. (2022) Mai, Q., Zhang, X., Pan, Y. and Deng, K. (2022).  A doubly enhanced em algorithm for model-based tensor clustering.  Journal of the American Statistical Association 117 2120–2134.
  * Mankad and Michailidis (2014) Mankad, S. and Michailidis, G. (2014).  Biclustering three-dimensional data arrays with plaid models.  Journal of Computational and Graphical Statistics 23 943–965.
  * Mitchell et al. (2008) Mitchell, T. M., Shinkareva, S. V., Carlson, A., Chang, K.-M., Malave, V. L., Mason, R. A. and Just, M. A. (2008).  Predicting human brain activity associated with the meanings of nouns.  Science 320 1191–1195.
  * Ning and Liu (2013) Ning, Y. and Liu, H. (2013).  High-dimensional semiparametric bigraphical models.  Biometrika 100 655–670.
  * OriGene (2020) OriGene (2020).  1500003o03rik (bc054733) mouse untagged clone.  https://www.origene.com/catalog/cdna-clones/expression-plasmids/mc206143/1500003o03rik-bc054733-mouse-untagged-clone, [Accessed: 2021-9-17].
  * Rand (1971) Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical association 66 846–850.
  * Sill et al. (2011) Sill, M., Kaiser, S., Benner, A. and Kopp-Schneider, A. (2011).  Robust biclustering by sparse singular value decomposition incorporating stability selection.  Bioinformatics 27 2089–2097.
  * Sun and Li (2019) Sun, W. W. and Li, L. (2019).  Dynamic tensor clustering.  Journal of the American Statistical Association 114 1894–1907.
  * Tan and Witten (2014) Tan, K. M. and Witten, D. M. (2014).  Sparse biclustering of transposable data.  Journal of Computational and Graphical Statistics 23 985–1008.
  * Vinh et al. (2009) Vinh, N. X., Epps, J. and Bailey, J. (2009).  Information theoretic measures for clusterings comparison: is a correction for chance necessary?  In Proceedings of the 26th annual international conference on machine learning.
  * Wang et al. (2019) Wang, M., Fischer, J. and Song, Y. S. (2019).  Three-way clustering of multi-tissue multi-individual gene expression data using semi-nonnegative tensor decomposition.  The annals of applied statistics 13 1103\.
  * Wang and Zeng (2019) Wang, M. and Zeng, Y. (2019).  Multiway clustering via tensor block models.  Advances in Neural Information Processing Systems 32 715–725.
  * Yin and Li (2012) Yin, J. and Li, H. (2012).  Model selection and estimation in the matrix normal graphical model.  Journal of multivariate analysis 107 119–140.
  * Zahn et al. (2007) Zahn, J. M., Poosala, S., Owen, A. B., Ingram, D. K., Lustig, A., Carter, A., Weeraratna, A. T., Taub, D. D., Gorospe, M. and Mazan-Mamczarz, K. (2007).  Agemap: a gene expression database for aging in mice.  PLoS genetics 3 e201.
  * Zhang et al. (2020) Zhang, Y., Levina, E. and Zhu, J. (2020).  Detecting overlapping communities in networks using spectral methods.  SIAM Journal on Mathematics of Data Science 2 265–283.
  * Zhao et al. (2016) Zhao, H., Wang, D. D., Chen, L., Liu, X. and Yan, H. (2016).  Identifying multi-dimensional co-clusters in tensors based on hyperplane detection in singular vector spaces.  PloS one 11 e0162293.

Supplementary Material

##  Appendix A Comparison with the Vectorized Approach in Bunea et al. (2020)

Model (1.1) can be viewed as the extension of the G-block model for clustering
a random vector in Bunea et al. (2020) to matrix valued data. If we vectorize
the matrix $X$, model (1.1) is equivalent to
$\textrm{vec}(X)=M\textrm{vec}(Z)+\textrm{vec}(\Gamma)$ with $M=B\otimes A$,
where $\textrm{vec}(X)$ denotes the vectorization of $X$, formed by stacking
the columns of $X$ into a single column vector, and $\otimes$ denotes the
Kronecker product. Thus, compared to Bunea et al. (2020) which allows $M$ to
be any unstructured $(pq)\times(K_{1}K_{2})$ membership matrix, we impose the
Kronecker product structure to the membership matrix $M$. While our model for
$\textrm{vec}(X)$ is more restrictive than Bunea et al. (2020), it actually
comes with two advantages for matrix clustering. First, as seen above, $A$ and
$B$ are interpreted as the membership matrices for the rows and columns.
Ignoring the Kronecker product structure and directly applying the model in
Bunea et al. (2020) would no longer produce interpretable results for matrix
clustering, as shown in Figure 3. Second, the Kronecker product of $A$ and $B$
provides a more parsimonious parametrization for the unknown membership matrix
$M$, leading to stronger theoretical guarantees on clustering.

Figure 3: $p=q=20$, $K_{1}=K_{2}=4$, and the row and column cluster sizes are
$(3,3,6,8)$ each. The data was generated with the matrix normal distribution
with a decaying Toeplitz matrix for the row and column covariance matrices and
a proportional noise variance setting. The left matrices are the cluster
results when using cord [Bunea et al. (2020)] on the vectorized matrix, and
the right matrices are the results when using our 2-STEP COD method. When
$n=50$, our method perfectly recovers the block cluster structure, but even
when the results are imperfect ($n=40$), the row and column clusters are still
easily interpretable, unlike the left results, which has cluster elements
scattered throughout the matrix without additional structure.

##  Appendix B Comparison of the Proposed Hierarchical Algorithm with the
Algorithm in Bunea et al. (2020)

Like shown in Figure 4, if C,D,E have not been clustered yet and A and B are
considered first, due to the nature of the algorithm in Bunea et al. (2020),
since the distance between C and A is less than $\alpha$, C will be clustered
with A and B even though C,D,E are closer to each other. If we implement a
hierarchical approach, the distances between all the pairs of points will be
considered at the same time, and a tree structure will naturally be
constructed. The two groups of {A,B} and {C,D,E} will be formed in the lower
part of the tree first, and even if the same $\alpha$ threshold value is used,
the two groups will not be merged together.

Figure 4: (a) Due to the nature of the algorithm from Bunea et al. (2020), A
and B are closest together and thus will be clustered before C,D,E are
considered. Then, since the distance between A and C is below the threshold
value $\alpha$, C will be incorrectly clustered with A and B. (b) With the
hierarchical algorithm, the distances between all the points are considered at
the same time and thus C will be correctly clustered with D and E.

##  Appendix C Conditions for $||X||_{W}=O(1)$

In this section, we establish the conditions for $||X||_{W}=O(1)$ for a
general class of weight matrices. For brevity, we will discuss the row
clustering case. Since our method is focused on using the columns to cluster
the rows, we will focus on the class of weight matrices that are generated by
arbitrary column membership matrices $\bar{B}$:

$\mathbf{W}_{(r)}(q):=\bigg{\\{}\bar{W}=\frac{1}{s}\bar{B}(\bar{B}^{T}\bar{B})^{-2}\bar{B}^{T}\bigg{|}\leavevmode\nobreak\
\bar{B}\text{ is a membership matrix}\in\mathbb{R}^{q\times s}\bigg{\\}},$

where $s$ is the number of clusters implied by $\bar{B}$. It can be shown that

$\displaystyle||X||_{\bar{W}}\leavevmode\nobreak\ $
$\displaystyle=\leavevmode\nobreak\
c\cdot\frac{\sqrt{K_{2}}}{s}\cdot\underset{1\leq a\leq
p}{\max}\sqrt{\Big{|}\Big{|}G^{T}G\Big{|}\Big{|}_{F}^{2}+\sum_{r=1}^{s}\Bigg{(}\frac{\sum_{j\in\bar{[r]}}\sigma^{2}_{aj}}{|\bar{[r]}|^{2}}\Bigg{)}^{2}}$
$\displaystyle\leq\leavevmode\nobreak\
c\cdot\frac{\sqrt{K_{2}}}{s}\cdot\sqrt{\Big{|}\Big{|}G^{T}G\Big{|}\Big{|}_{F}^{2}+s\cdot\Bigg{(}\frac{\sigma^{2}_{\max}}{\bar{m}_{q}}\Bigg{)}^{2}},$

where $G=B^{T}\bar{B}(\bar{B}^{T}\bar{B})^{-1}$ and $\bar{m}_{q}$ denotes the
smallest column cluster size implied by $\bar{B}$. So,

$||X||_{\bar{W}}\leavevmode\nobreak\ \lesssim\leavevmode\nobreak\ O(1)\hskip
14.22636pt\text{ if both }\leavevmode\nobreak\ \leavevmode\nobreak\
\begin{cases}\leavevmode\nobreak\ \leavevmode\nobreak\
||G^{T}G||_{F}^{2}\leavevmode\nobreak\ &\lesssim\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \frac{s^{2}}{K_{2}}\\\
\leavevmode\nobreak\ \leavevmode\nobreak\
\sigma^{2}_{\max}\leavevmode\nobreak\ &\lesssim\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\bar{m}_{q}\cdot\sqrt{\frac{s}{K_{2}}}\end{cases}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \text{hold}.$

Note that when $\bar{B}=B$, then $G=I_{s}=I_{K_{2}}$ and the condition for
$||G^{T}G||_{F}^{2}$ is automatically satisfied. We expect that the condition
is also satisfied for small perturbations $\bar{B}\approx B$, where $s\asymp
K_{2}$.

##  Appendix D Derivation for Remark 3.2

We consider the setting where the column cluster has equal size (i.e.,
$K_{2}M_{q}=q$) and $\sigma_{ab}^{2}=\sigma^{2}$. First, we have shown in
Section 2.2 that the off diagonal entries of $\Sigma_{p,W_{I}}$ and
$\Sigma_{p,W_{O}}$ are the same, that is
$(\Sigma_{p,W_{I}})_{jk}=(\Sigma_{p,W_{O}})_{jk},\leavevmode\nobreak\ j\neq
k$. Then we have
$\mathrm{MCOD}(\Sigma_{p,W_{I}})=\mathrm{MCOD}(\Sigma_{p,W_{O}}):=\mathrm{MCOD}$.
For $W_{I}$, we can show that

$\displaystyle\frac{\mathrm{MCOD}(\Sigma_{p,W_{I}})}{\|X\|_{W_{I}}}$
$\displaystyle\geq\frac{\mathrm{MCOD}}{\sqrt{\frac{1}{K_{2}}}\max_{k}||\operatorname{{\rm
Var}}(Z_{k\cdot})||_{F}+\sqrt{\frac{K_{2}}{q}}\sigma^{2}}$ $\displaystyle\geq
C\min\bigg{(}\frac{\mathrm{MCOD}}{\sqrt{\frac{1}{K_{2}}}\max_{k}||\operatorname{{\rm
Var}}(Z_{k\cdot})||_{F}},\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{\mathrm{MCOD}}{\sqrt{\frac{K_{2}}{q}}\sigma^{2}}\bigg{)},$

where $C$ is a constant. Hence, the cluster separation condition (A1) for
$W_{I}$ is implied by

$\frac{\mathrm{MCOD}}{\sqrt{\frac{1}{K_{2}}}\max_{k}||\operatorname{{\rm
Var}}(Z_{k\cdot})||_{F}}\geq c_{0}/C\eta,\leavevmode\nobreak\
\leavevmode\nobreak\
\frac{\mathrm{MCOD}}{\sqrt{\frac{K_{2}}{q}}\sigma^{2}}\geq c_{0}/C\eta.$

Similarly, the cluster separation condition (A1) for $W_{O}$ is implied by

$\frac{\mathrm{MCOD}}{\sqrt{\frac{1}{K_{2}}}\max_{k}||\operatorname{{\rm
Var}}(Z_{k\cdot})||_{F}}\geq c_{0}/C\eta,\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{\mathrm{MCOD}}{{\frac{K_{2}}{q}}\sigma^{2}}\geq
c_{0}/C\eta.$

In terms of the noise level $\sigma^{2}$, the condition for $W_{I}$ is
$\sigma^{2}=O(\frac{\mathrm{MCOD}\sqrt{q}}{\eta\sqrt{K_{2}}})$, whereas the
condition for $W_{O}$ is $\sigma^{2}=O(\frac{\mathrm{MCOD}\cdot q}{\eta
K_{2}})$. The latter is clearly much weaker. Thus, when $q\gg K_{2}$, our
algorithm with the optimal weight $W_{O}$ attains clustering consistency in
the presence of a larger noise level $\sigma^{2}$. This is the benefit of
using the optimal weight in our algorithm.

##  Appendix E Minimax Lower Bound with a Perturbed $\bar{B}$

Here we present a lower bound result when the cluster separation metric in the
definition of the parameter space is defined based on an arbitrary column
membership matrix $\bar{B}$. For example, the special cases $\bar{B}=I_{q}$
and $\bar{B}=B$ are considered.

Assume that we have an arbitrary column membership matrix $\bar{B}$, which is
fixed and can be different from the true membership matrix $B$. Let
$\bar{W}=\bar{B}(\bar{B}^{T}\bar{B})^{-2}\bar{B}^{T}/s$ where $s$ denotes the
number of clusters implied by $\bar{B}$. We define the following parameter
space

$M_{\bar{W}}=\big{\\{}\Sigma\in\mathbb{R}^{pq\times pq}|\leavevmode\nobreak\
X\leavevmode\nobreak\ \textrm{satisfies model}\leavevmode\nobreak\
(1.1),\leavevmode\nobreak\ \leavevmode\nobreak\
\text{MCOD}(\Sigma_{p,\bar{W}})/\|X\|_{\bar{W}}\geq\eta\big{\\}},$

We now present a general lower bound for clustering over the parameter space
$M_{\bar{W}}$.

######  Theorem E.1.

(Minimax Lower Bound with a Perturbed $\bar{B}$)
For $K_{1}\geq 3$, there exists a positive constant $c$ such that for any
$\eta$

$\displaystyle 0\leavevmode\nobreak\ \leavevmode\nobreak\
\leq\leavevmode\nobreak\ \leavevmode\nobreak\ \eta\leavevmode\nobreak\
\leavevmode\nobreak\ <\leavevmode\nobreak\ \leavevmode\nobreak\
c\sqrt{\frac{\log
p}{nK_{2}}}\cdot\frac{||G||_{F}^{2}}{\sqrt{K_{2}}\cdot\sqrt{\big{|}\big{|}G^{T}G\big{|}\big{|}_{F}^{2}+\Big{(}\frac{q}{K_{2}}\Big{)}^{2}\cdot\frac{s}{\bar{m}_{q}^{2}}}}$
(E.1)

we have

$\displaystyle\underset{\widehat{\mathcal{G}}}{\inf}\leavevmode\nobreak\
\underset{\Sigma\in M_{\bar{W}}}{\sup}\leavevmode\nobreak\
\mathbb{P}_{\Sigma}(\widehat{\mathcal{G}}\neq\mathcal{G})\leavevmode\nobreak\
\geq\leavevmode\nobreak\ \frac{1}{7},$

where $\bar{m}_{q}$ denotes the smallest column cluster size implied by
$\bar{B}$ and $G=B^{T}\bar{B}(\bar{B}^{T}\bar{B})^{-1}$. The infimum is taken
over all possible estimators of $\mathcal{G}$.

This theorem shows that, if we define MCOD based on a perturbed $\bar{B}$, the
rate of the cluster separation (E.1) depends on $G$ and $\bar{m}_{q}$. If
$\bar{B}$ is near perfect, $G\approx I_{K_{2}}$, $G^{T}G\approx I_{K_{2}}$,
$s\approx K_{2}$ and $\bar{m}_{q}\approx m_{q}=\frac{q}{K_{2}}$ in the
construction of the lower bound. Then, $||G||_{F}^{2}$, $||G^{T}G||_{F}^{2}$
and $(\frac{q}{K_{2}})^{2}\cdot\frac{s}{\bar{m}_{q}^{2}}$ will all be close to
order $K_{2}$ and the bound in (E.1) becomes $\sqrt{\frac{\log p}{nK_{2}}}$.
In contrast, if we use $W_{I}=\frac{1}{q}I_{q}$, then $G$ is a $K_{2}\times q$
matrix with row vectors $1^{T}_{m_{q}^{(i)}}$ as blocks on the diagonal. This
implies that $G^{T}G$ is a $q\times q$ matrix with
$1_{m_{q}^{(i)}}1^{T}_{m_{q}^{(i)}}$ as square blocks in the diagonal,
$||G||_{F}^{2}=q$ and
$||G^{T}G||_{F}^{2}=K_{2}\cdot(\frac{q}{K_{2}})^{2}=\frac{q^{2}}{K_{2}}$ in
our construction. The lower bound becomes $\sqrt{\frac{\log p}{nq}}$. These
two cases show how the lower bound depends on the imperfect column cluster
structure $\bar{B}$. The proof of this lower bound can be found in Section
H.2.2 of the Supplementary Material.

##  Appendix F Supplementary Material for Section 5

###  F.1 Theoretical Guarantees for Algorithm 2

####  F.1.1 Consistency

We present the consistency theorem for the iterative one-step Algorithm 2 that
clusters both the rows and the columns.

######  Theorem F.1.

(Consistency with One-step Hierarchical COD with the Optimal Weight)
Under the model $X=AZB^{T}+\Gamma$, assume that $\textrm{vec}(X)$ is
multivariate Gaussian, $\log p=o(n)$, $\log q=o(n)$ and the following
conditions hold:

  * (R)

Row Separation Condition:

$\text{MCOD}(\Sigma_{p,W_{O,(r)}})/\|X\|_{W_{O,(r)}}\leavevmode\nobreak\
\leavevmode\nobreak\ >\leavevmode\nobreak\ \leavevmode\nobreak\
c_{0,(r)}\cdot\eta_{(r)}$

where $c_{0,(r)}\geq 4$ is an arbitrary constant and $\leavevmode\nobreak\
\eta_{(r)}\geq c_{1,(r)}\sqrt{\frac{\log p}{nK_{2}}}\leavevmode\nobreak\ $ for
a universal constant $c_{1,(r)}>0$.

  * (C)

Column Separation Condition:

$\text{MCOD}(\Sigma_{q,W_{O,(c)}})/\|X\|_{W_{O,(c)}}\leavevmode\nobreak\
\leavevmode\nobreak\ >\leavevmode\nobreak\ \leavevmode\nobreak\
c_{0,(c)}\cdot\eta_{(c)}$

where $c_{0,(c)}\geq 4$ is an arbitrary constant and $\leavevmode\nobreak\
\eta_{(c)}\geq c_{1,(c)}\sqrt{\frac{\log q}{nK_{1}}}\leavevmode\nobreak\ $ for
a universal constant $c_{1,(c)}>0$.

  * (S)

Stability Condition:

$\bigg{[}\Big{\\{}\text{MCOD}(\Sigma_{p,\widehat{W}_{O,(r)}})\Big{\\}}\Big{/}\|X\|_{\widehat{W}_{O,(r)}}\bigg{]}\Bigg{/}\bigg{[}\Big{\\{}\text{MCOD}(\Sigma_{p,W_{O,(r)}})\Big{\\}}\Big{/}\|X\|_{W_{O,(r)}}\bigg{]}>\frac{4}{c_{0,(r)}}$

where $c_{0,(r)}$ is defined in condition (R).

Then using our Algorithm 2 with the threshold
$\alpha_{(r)}=2\eta_{(r)}\cdot\|X\|_{\widehat{W}_{O,(r)}}$ for the row
clustering step and the threshold
$\alpha_{(c)}=2\eta_{(c)}\cdot\|X\|_{\widehat{W}_{O,(c)}}$ for the column
clustering step, we obtain perfect cluster recovery (i.e.,
$\widehat{\mathcal{G}}^{(r)}=\mathcal{G}^{(r)}$ and
$\widehat{\mathcal{G}}^{(c)}=\mathcal{G}^{(c)}$) with probability greater than
$1-\frac{c_{3}}{p}$ for some constant $c_{3}>0$.

The row and column separation conditions (R) and (C) are similar to the
cluster separation condition (A1) in Theorem 3.1. The stability condition (S)
is only imposed for row clustering, which can be viewed as the condition for
the initial value $\bar{B}$ in Step 1 of Algorithm 2. As shown in Proposition
4.2, the initial value $\bar{B}=I_{q}$ satisfies condition (S) under the
matrix normal model with mild additional assumptions. Indeed, by Theorem 3.1,
conditions (S) and (R) imply perfect row cluster recovery (i.e.,
$\widehat{A}=A$) with high probability, which further implies that
$\widehat{W}_{O,(c)}=W_{O,(c)}$. Thus, the stability condition for column
clustering automatically holds, which in turn guarantees perfect column
cluster recovery with high probability. The proof is very similar to the proof
for Theorem 3.1 and will be omitted for brevity.

####  F.1.2 Minimax Optimality

In the following, we establish a lower bound result for both row and column
clustering which naturally incorporates the uncertainty in estimating the
membership matrices $A$ and $B$ simultaneously. Following the notation used in
Theorem F.1, we define parameter spaces

$M_{O,(r)}=\big{\\{}\Sigma\in\mathbb{R}^{pq\times pq}|\leavevmode\nobreak\
X\leavevmode\nobreak\ \textrm{satisfies model}\leavevmode\nobreak\
(1.1),\leavevmode\nobreak\ \leavevmode\nobreak\
\text{MCOD}(\Sigma_{p,W_{O,(r)}})/\|X\|_{W_{O,(r)}}\geq\eta_{(r)}\big{\\}},$

and

$M_{O,(c)}=\big{\\{}\Sigma\in\mathbb{R}^{pq\times pq}|\leavevmode\nobreak\
X\leavevmode\nobreak\ \textrm{satisfies model}\leavevmode\nobreak\
(1.1),\leavevmode\nobreak\ \leavevmode\nobreak\
\text{MCOD}(\Sigma_{q,W_{O,(c)}})/\|X\|_{W_{O,(c)}}\geq\eta_{(c)}\big{\\}}$

for the rows and columns, respectively. To study the minimax lower bound for
both row and column clustering, we consider the parameter space $M_{O,(r)}\cap
M_{O,(c)}$. The following theorem provides the lower bound for clustering over
this parameter space $M_{O,(r)}\cap M_{O,(c)}$.

######  Theorem F.2.

For $K_{1},K_{2}\geq 3$, there exists a positive constant $c$ such that, for
any $\eta_{(r)}$ and $\eta_{(c)}$

$\displaystyle 0\leq\eta_{(r)}<c\sqrt{\frac{\log
p}{nK_{2}}}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\textrm{or}\leavevmode\nobreak\ \leavevmode\nobreak\
0\leq\eta_{(c)}<c\sqrt{\frac{\log q}{nK_{1}}},$

we have

$\displaystyle\underset{(\widehat{\mathcal{G}}^{(r)},\widehat{\mathcal{G}}^{(c)})}{\inf}\leavevmode\nobreak\
\leavevmode\nobreak\ \underset{\Sigma\in M_{O,(r)}\cap
M_{O,(c)}}{\sup}\leavevmode\nobreak\ \leavevmode\nobreak\
\mathbb{P}_{\Sigma}\Big{(}(\widehat{\mathcal{G}}^{(r)},\widehat{\mathcal{G}}^{(c)})\neq(\mathcal{G}^{(r)},\mathcal{G}^{(c)})\Big{)}\leavevmode\nobreak\
\leavevmode\nobreak\ \geq\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{1}{7},$

where the infimum is taken over all possible estimators of
$(\mathcal{G}^{(r)},\mathcal{G}^{(c)})$.

This theorem shows that we need $\eta_{(r)}\geq c\sqrt{\frac{\log p}{nK_{2}}}$
and $\eta_{(c)}\geq c\sqrt{\frac{\log q}{nK_{1}}}$ to attain both perfect row
and column clustering. Together with Theorem F.1, we obtain that, as long as
the initial estimate $\bar{B}$ satisfies the stability condition (i.e., it
falls into a contraction region), the One-Step Hierarchical Algorithm with the
Optimal Weight (Algorithm 2) is minimax optimal for row and column clustering.
Recall that this has been verified for the matrix normal model with
$\bar{B}=I_{q}$ in Proposition 4.2. Finally, we note that since the above
lower bound is concerned with both row and column clustering, the uncertainty
in estimating both membership matrices $A$ and $B$ is taken into account. The
proof of this theorem is in Section H.2.3 of the Supplementary Material.

###  F.2 Two-Step Hierarchical Algorithm with the Optimal Weight

The two-step algorithm is shown in Algorithm 3. In this algorithm, we ignore
the sample splitting step for simplicity.

  1. (0a)

Set the initial value $\bar{B}=I_{q}$.

  2. (0b)

Apply Algorithm 1 with
$\widehat{\Sigma}_{p,\widehat{W}_{O_{1},(r)}}=\frac{1}{|D_{2}|}\sum_{i\in
D_{2}}X^{(i)}\widehat{W}_{O_{1},(r)}X^{(i)T}$ to cluster the

rows of $X$, where
$\widehat{W}_{O_{1},(r)}=\bar{B}(\bar{B}^{T}\bar{B})^{-2}\bar{B}^{T}/s$ and
$s$ denotes the estimated number

of clusters from $\bar{B}$. Obtain the resulting row cluster
$\widehat{\mathcal{G}}^{(r)}$ or equivalently the membership matrix
$\widehat{A}_{1}$.

* (1a)

Compute the estimate of the optimal column weight
$\widehat{W}_{O,(c)}=\widehat{A}_{1}(\widehat{A}_{1}^{T}\widehat{A}_{1})^{-2}\widehat{A}_{1}^{T}/t$
where $t$ is the estimated number of row clusters in $\widehat{A}_{1}$. (1b)
Apply Algorithm 1 with
$\widehat{\Sigma}_{q,\widehat{W}_{O,(c)}}=\frac{1}{n}\sum_{i=1}^{n}X^{(i)T}\widehat{W}_{O,(c)}X^{(i)}$
to cluster the columns

of $X$ and find the estimator $\widehat{B}_{1}$.

* (2a)

Compute the estimate of the optimal row weight
$\widehat{W}_{O_{2},(r)}=\widehat{B}_{1}(\widehat{B}_{1}^{T}\widehat{B}_{1})^{-2}\widehat{B}_{1}^{T}/s_{1}$,

where $s_{1}$ denotes the estimated number of column clusters in
$\widehat{B}_{1}$.

* (2b)

Apply Algorithm 1 with
$\widehat{\Sigma}_{p,\widehat{W}_{O_{2},(r)}}=\frac{1}{n}\sum_{i=1}^{n}X^{(i)}\widehat{W}_{O_{2},(r)}X^{(i)T}$
to cluster the rows

of $X$. Obtain the final estimator $\widehat{A}_{2}$.

Algorithm 3 Two-step Hierarchical Algorithm with the Optimal Weight

###  F.3 Simulation Results With and Without Sample Splitting

We present the results for an additional simulation study that highlights the
difference of our method with and without sample splitting in Table 1. We have
$n$ i.i.d. copies of $30\times 30$ matrices with $K_{1}=K_{2}=4$ and
moderately unbalanced cluster sizes of $4,6,9,11$ for both the rows and
columns. The decay rate for the Toeplitz covariance matrices is -0.2 and 0.2
for the rows and columns, respectively. We consider the proportional noise
variance setting from the main paper. We vary $n$ from $20$ up to $100$ in
increments of $20$ and the ARI values for 2-STEP COD are recorded. We see
that, for $n$ relatively small (say $n\leq 60$), our 2-STEP COD without data
splitting performs significantly better than the method using data splitting.
In addition, both methods yield very high clustering accuracy when $n$ is
large enough (say $n\geq 100$). Thus, in practice, we recommend using our
2-STEP COD without data splitting, especially when the sample size is small or
moderate.

|  | Data Split (Row) | No Data Split (Row) | Data Split (Col) | No Data Split (Col)
---|---|---|---|---|---
$n$ | 20 | 0 | 0.4984 | 0 | 0.2723
40 | 0.0639 | 0.9939 | 0.0712 | 0.9562
| 60 | 0.5642 | 1 | 0.1647 | 0.9979
| 80 | 0.9685 | 1 | 0.6834 | 0.9934
| 100 | 0.9849 | 1 | 0.9528 | 0.9962

Table 1: The ARI values obtained from the 2-STEP COD Algorithm with $p=q=30$,
$K_{1}=K_{2}=4$, moderately unbalanced cluster sizes with proportional noise
variance under the $n=20,40,60,80,100$ setting with and without data
splitting.

###  F.4 Data-Driven Tuning Parameter Selection Process for $\alpha$

In the following, we describe a data-driven selection method for $\alpha$ that
was also used in Bunea et al. (2020). To fix the notation, we consider our
Algorithm 1 with some sample covariance matrix $\widehat{\Sigma}_{p}$. We
assume the data are standardized. The method can be similarly applied to the
one-step, two-step algorithm with the optimal weight and the “naive”
algorithm. The steps are outlined in Algorithm 6.

We first split the data into two, $D_{1}$, $D_{2}$, and calculate
$\widehat{\Sigma}_{p}^{(1)}$, $\widehat{\Sigma}_{p}^{(2)}$, respectively. For
each tuning parameter $\alpha_{l}$ in the grid, we perform our algorithm on
$\widehat{\Sigma}_{p}^{(1)}$ to get a cluster structure $\mathcal{G}_{l}$. We
then take the average of all the non-diagonal elements in the cluster blocks
of $\widehat{\Sigma}_{p}^{(1)}$ via the smoothing operator
$\Upsilon(\widehat{\Sigma}_{p}^{(1)},\mathcal{G}_{l})$. Finally we calculate
the Frobenius loss of $\Upsilon(\Sigma_{p}^{(1)},\mathcal{G}_{l})$ and
$\Sigma_{p}^{(2)}$, and choose $\alpha_{l}$ that yields the smallest value.

  1. 1.

Split the data into two: $D_{1}$ and $D_{2}$

  2. 2.

Using $D_{1}$, calculate $\widehat{\Sigma}_{p}^{(1)}$.

  3. 3.

Using $D_{2}$, calculate $\widehat{\Sigma}_{p}^{(2)}$.

  4. 4.

For $r>1$ and each value $\alpha_{l}$ on a grid $(l=1,...,r)$, perform
Algorithm 1 with $\widehat{\Sigma}_{p}^{(1)}$

to get a row cluster structure $\mathcal{G}_{l}^{(r)}$.

* 5.

Perform the smoothing operator
$\Upsilon(\widehat{\Sigma}_{p}^{(1)},\mathcal{G}_{l}^{(r)})$ where $\Upsilon$
is defined as the following:

$[\Upsilon(\widehat{\Sigma}_{p},\mathcal{G}^{(r)})]_{ab}=\begin{cases}\Big{|}G_{k}^{(r)}\Big{|}^{-1}\bigg{(}\Big{|}G_{k}^{(r)}\Big{|}-1\bigg{)}^{-1}\underset{i,j\in
G_{k},\leavevmode\nobreak\ i\neq
j}{\sum}\big{[}\widehat{\Sigma}_{p}\big{]}_{ij}&\text{if}\leavevmode\nobreak\
a\neq b,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ k=k^{\prime}\\\
\Big{|}G_{k}^{(r)}\Big{|}^{-1}\Big{|}G_{k^{\prime}}^{(r)}\Big{|}^{-1}\underset{i\in
G_{k},\leavevmode\nobreak\ j\in
G_{k^{\prime}}}{\sum}\big{[}\widehat{\Sigma}_{p}\big{]}_{ij}&\text{if}\leavevmode\nobreak\
a\neq b,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ k\neq
k^{\prime}\\\ 1&\text{if}\leavevmode\nobreak\ a=b.\end{cases}$

* 6.

Our data dependent tuning parameter for the threshold is:

$\widehat{\alpha}=\underset{\alpha_{l}}{\mathrm{argmin}}\leavevmode\nobreak\
L\bigg{(}\Upsilon(\widehat{\Sigma}_{p}^{(1)},\mathcal{G}_{l}^{(r)}),\leavevmode\nobreak\
\widehat{\Sigma}_{p}^{(2)}\bigg{)},$

where $L(A,B):=||A-B||_{F}$ is the Frobenius loss.

Algorithm 4 A Data-Driven Tuning Parameter Selection Process

##  Appendix G Supplementary Material for Section 6

###  G.1 The Dependent Noise Model

####  G.1.1 Discussion of Theorem 6.1

Under model (1.1) with dependent noise elements, since
$\Sigma_{p,W}=\mathbb{E}(AZB^{T}WBZ^{T}A^{T})+\mathbb{E}(\Gamma W\Gamma^{T})$,
if $\leavevmode\nobreak\ a\underset{\mathcal{G}}{\sim}b$,

$\displaystyle\text{COD}_{\Sigma_{p,W}}(a,b)\leavevmode\nobreak\
\leavevmode\nobreak\ $ $\displaystyle=\leavevmode\nobreak\
\leavevmode\nobreak\ \underset{c\neq a,b}{\max}\leavevmode\nobreak\
\Big{|}\big{[}\Sigma_{p,W}\big{]}_{ac}-\big{[}\Sigma_{p,W}\big{]}_{bc}\Big{|}$
$\displaystyle=\leavevmode\nobreak\ \leavevmode\nobreak\ \underset{c\neq
a,b}{\max}\leavevmode\nobreak\
\Big{|}\big{[}\mathbb{E}(AZB^{T}WBZ^{T}A^{T})\big{]}_{ac}-\big{[}\mathbb{E}(AZB^{T}WBZ^{T}A^{T})\big{]}_{bc}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \hskip
156.49014pt+\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{ac}-\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{bc}\Big{|}$ $\displaystyle=\leavevmode\nobreak\
\leavevmode\nobreak\ \underset{c\neq a,b}{\max}\leavevmode\nobreak\
\Big{|}\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{ac}-\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{bc}\Big{|}$ $\displaystyle\leq\leavevmode\nobreak\
\leavevmode\nobreak\ \gamma(\Sigma_{p,W})$

and if $\leavevmode\nobreak\ a\underset{\mathcal{G}}{\not\sim}b$,

$\displaystyle\text{COD}_{\Sigma_{p,W}}(a,b)\leavevmode\nobreak\
\leavevmode\nobreak\ $ $\displaystyle=\leavevmode\nobreak\
\leavevmode\nobreak\ \underset{c\neq a,b}{\max}\leavevmode\nobreak\
\Big{|}\big{[}\mathbb{E}(AZB^{T}WBZ^{T}A^{T})\big{]}_{ac}-\big{[}\mathbb{E}(AZB^{T}WBZ^{T}A^{T})\big{]}_{bc}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \hskip
156.49014pt+\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{ac}-\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{bc}\Big{|}$ $\displaystyle\geq\leavevmode\nobreak\
\leavevmode\nobreak\ \underset{c\neq a,b}{\max}\leavevmode\nobreak\
\Big{|}\big{[}\mathbb{E}(AZB^{T}WBZ^{T}A^{T})\big{]}_{ac}-\big{[}\mathbb{E}(AZB^{T}WBZ^{T}A^{T})\big{]}_{bc}\Big{|}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \hskip
156.49014pt-\underset{c\neq a,b}{\max}\leavevmode\nobreak\
\Big{|}\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{ac}-\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{bc}\Big{|}$ $\displaystyle\geq\leavevmode\nobreak\
\leavevmode\nobreak\ \text{MCOD*}(\Sigma_{p,W})\leavevmode\nobreak\
\leavevmode\nobreak\ -\leavevmode\nobreak\ \leavevmode\nobreak\
\underset{c\neq a,b}{\max}\leavevmode\nobreak\ \Big{|}\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{ac}-\big{[}\mathbb{E}(\Gamma
W\Gamma^{T})\big{]}_{bc}\Big{|}$ $\displaystyle\geq\leavevmode\nobreak\
\leavevmode\nobreak\ \text{MCOD*}(\Sigma_{p,W})\leavevmode\nobreak\
\leavevmode\nobreak\ -\leavevmode\nobreak\ \leavevmode\nobreak\
\gamma(\Sigma_{p,W}),$

where
$\text{MCOD*}(\Sigma_{p,W})=\underset{a\underset{\mathcal{G}}{\not\sim}b}{\min}\leavevmode\nobreak\
\underset{c\neq a,b}{\max}\leavevmode\nobreak\
\Big{|}\big{[}\mathbb{E}(AZB^{T}WBZ^{T}A^{T})\big{]}_{ac}-\big{[}\mathbb{E}(AZB^{T}WBZ^{T}A^{T})\big{]}_{bc}\Big{|}$.
In other words, MCOD* is the counterpart to MCOD that only considers the
signal (the first component) from
$\Sigma_{p,W}=\mathbb{E}(AZB^{T}WBZ^{T}A^{T})+\mathbb{E}(\Gamma W\Gamma^{T})$.

Comparing the population quantity $\text{COD}_{\Sigma_{p,W}}(a,b)$ in the two
cases $\leavevmode\nobreak\ a\underset{\mathcal{G}}{\sim}b$ and
$\leavevmode\nobreak\ a\underset{\mathcal{G}}{\not\sim}b$, it is apparent that
if the signal term is strong enough in the following sense

$\text{MCOD*}(\Sigma_{p,W})\leavevmode\nobreak\ \leavevmode\nobreak\
>\leavevmode\nobreak\ \leavevmode\nobreak\ 2\cdot\gamma(\Sigma_{p,W}),$

then even when the noise variables are dependent, the COD measures are well
separated and on the population level, the COD method is still viable for
clustering.

The proof of this extended model can be constructed using the above logic and
the proof from Theorem 3.1 which is in Section H.1.1 of the Supplementary
Material. The remaining parts will be omitted for brevity.

####  G.1.2 A Slightly Different Measure, $\gamma_{s}(\Sigma_{p,W})$

We can swap out the definition of $\gamma(\Sigma_{p,W})$ in (6.1) for a
slightly more restrictive but more intuitive defintion:

$\displaystyle\gamma_{s}(\Sigma_{p,W})\leavevmode\nobreak\
:=\leavevmode\nobreak\ 2\underset{1\leq a\neq b\leq
p}{\max}\leavevmode\nobreak\ \leavevmode\nobreak\
\Big{|}\big{[}\mathbb{E}(\Gamma W\Gamma^{T})\big{]}_{ab}\Big{|}$ (G.1)

This $\gamma_{s}(\Sigma_{p,W})$ can be used in place of the
$\gamma(\Sigma_{p,W})$ in the above theorem since the former is an upper bound
for the latter. $\gamma_{s}(\Sigma_{p,W})$ gives information on the largest
off-diagonal entry in $\mathbb{E}(\Gamma W\Gamma^{T})$, i.e. the largest
weighted covariance between two different noise variables. The relationship
between this weighted covariance and the actual covariance between noise
variables depends on the weight. For example, with the optimal weight $W_{O}$,
$\gamma_{s}(\Sigma_{p,W_{O}})$ can be upper bounded with
$2\cdot\sigma_{\text{offmax}}^{2}$, where
$\sigma_{\text{offmax}}^{2}:=\underset{(i,j)\neq(i^{\prime},j^{\prime})}{\max}\text{Cov}(\Gamma_{ij},\Gamma_{i^{\prime}j^{\prime}})$,
the maximum of the unweighted covariance between two different noise
variables.

###  G.2 Nested Clustering to Incorporate Mean Information

####  G.2.1 The Generalized Model

To generalize the proposed method to account for both mean and covariance
information in clustering, we extend our latent variable model to

$\displaystyle X=M+AZB^{T}+\Gamma,$ (G.2)

where $M=\mathbb{E}(X)$, $Z$ and $\Gamma$ are mean 0 random matrices. Like
before, $AZB^{T}$ induces the row and column clustering structures based on
the covariance of $X$. Specifically, $A\in\mathbb{R}^{p\times K_{1}}$ and
$B\in\mathbb{R}^{q\times K_{2}}$ are the unknown membership matrices for the
rows and columns, respectively. We define the row clusters as

$\mathcal{G}^{(r),[2]}=\\{G_{1}^{(r),[2]},...,G_{K_{1}}^{(r),[2]}\\},\leavevmode\nobreak\
\leavevmode\nobreak\ \textrm{where}\leavevmode\nobreak\ \leavevmode\nobreak\
G_{k}^{(r),[2]}=\\{a:A_{ak}=1\\}$ (G.3)

for any $1\leq k\leq K_{1}$. The column clusters $\mathcal{G}^{(c),[2]}$ can
be defined similarly. To incorporate the mean information, we further assume
that the matrix $M$ induces the row and column clustering structures based on
the mean of $X$. In particular, we assume $M=\widetilde{A}T\widetilde{B}^{T}$,
where $\widetilde{A}\in\mathbb{R}^{p\times\widetilde{K}_{1}}$ and
$B\in\mathbb{R}^{q\times\widetilde{K}_{2}}$ are the unknown membership
matrices for the rows and columns based on the mean information, and
$T\in\mathbb{R}^{\widetilde{K}_{1}\times\widetilde{K}_{2}}$ is an unknown
deterministic matrix. The corresponding row clusters are defined as

$\mathcal{G}^{(r),[1]}=\\{G_{1}^{(r),[1]},...,G_{\widetilde{K}_{1}}^{(r),[1]}\\},\leavevmode\nobreak\
\leavevmode\nobreak\ \textrm{where}\leavevmode\nobreak\ \leavevmode\nobreak\
G_{k}^{(r),[1]}=\\{a:\widetilde{A}_{ak}=1\\}$ (G.4)

for any $1\leq k\leq\widetilde{K}_{1}$. The column clusters
$\mathcal{G}^{(c),[1]}$ can be defined similarly. To link these two row
clusters $\mathcal{G}^{(r),[2]}$ and $\mathcal{G}^{(r),[1]}$, we assume
$\mathcal{G}^{(r),[2]}$ is nested inside $\mathcal{G}^{(r),[1]}$.

Definition 1. (Li et al., 2010) A clustering $\mathcal{G}^{[2]}$ with
$K^{\prime}$ clusters is said to be nested inside another clustering
$\mathcal{G}^{[1]}$ with $K$ clusters if:

  1. 1.

(Hierarchical Structure) For any cluster
$G^{[2]}_{j}\in\mathcal{G}^{[2]},(1\leq j\leq K^{\prime})$, there is a cluster
$G^{[1]}_{i}\in\mathcal{G}^{[1]},(1\leq i\leq K)$ such that
$G^{[2]}_{j}\subseteq G^{[1]}_{i}$.

  2. 2.

(Proper Subset Structure) There exists at least one cluster in
$\mathcal{G}^{[2]}$, (i.e. $G^{[2]}_{j*}$), which satisfies
$G^{[2]}_{j*}\subset G^{[1]}_{i*}$ and $G^{[2]}_{j*}\neq G^{[1]}_{i*}$, for
some cluster $G^{[1]}_{i*}\in\mathcal{G}^{[1]}$.

In other words, the clustering $\mathcal{G}^{(r),[2]}$ based on the covariance
information provides a more refined partition on top of the initial clustering
$\mathcal{G}^{(r),[1]}$ from the mean information. For the ease of
interpretation, we only consider the case that $\mathcal{G}^{(r),[2]}$ is
nested in $\mathcal{G}^{(r),[1]}$. This assumption can be relaxed or even
removed by defining clusters at different levels (e.g., clustering from the
mean and covariance) and new ways of combining the clusters.

####  G.2.2 Nested Clustering (Algorithm 5)

To recover the nested cluster structure, we propose a two-step nested
clustering algorithm, in which a mean-based clustering method is implemented
first, and then on each cluster, our covariance-based method is applied to
capture the finer, more intricate relationships within each broad cluster.

There are many existing mean-based clustering methods in the literature. Here
we use SparseBC from Tan and Witten (2014) as the mean-based clustering
method. The resulting clusterings are denoted by
$\widehat{\mathcal{G}}^{(r),[1]}$ and $\widehat{\mathcal{G}}^{(c),[1]}$. Then,
our proposed two-step Algorithm 3 is implemented to get the second layer
cluster structures within each first layer cluster. More specifically, for row
clustering, we apply the two-step Algorithm 3 to each
$\big{|}\widehat{G}^{(r),[1]}_{j}\big{|}\times q$ submatrix corresponding to
the variables in $\widehat{G}^{(r),[1]}_{j}$ for $1\leq
j\leq\big{|}\widehat{\mathcal{G}}^{(r),[1]}\big{|}$. Similarly, for column
clustering, we apply our two-step Algorithm 3 to each
$p\times\big{|}\widehat{G}^{(c),[1]}_{j}\big{|}$ submatrix corresponding to
the variables in $\widehat{G}^{(c),[1]}_{j}$ for $1\leq
j\leq\big{|}\widehat{\mathcal{G}}^{(c),[1]}\big{|}$. In this way, we can use
mean information to reduce the dimension of the matrices, to which our two-
step Algorithm 3 is applied. Thus, both the mean and the covariance
information is utilized in deriving the final clustering result. We summarize
this method in Algorithm 5.

  1. (1)

Perform mean-based clustering (SparseBC) to get the first layer cluster
structures $\widehat{\mathcal{G}}^{(r),[1]}$ and
$\widehat{\mathcal{G}}^{(c),[1]}$

  2. (2a)

Perform Algorithm 3 on each $\Big{|}\widehat{G}_{j}^{(r),[1]}\Big{|}\times q$
submatrix $(1\leq j\leq\Big{|}\widehat{\mathcal{G}}^{(r),[1]}\Big{|})$

  3. (2b)

Perform (the col. version of) Algorithm 3 on each
$p\times\Big{|}\widehat{G}_{j}^{(c),[1]}\Big{|}$ submatrix $(1\leq
j\leq\Big{|}\widehat{\mathcal{G}}^{(c),[1]}\Big{|})$

  4. (3a)

Combine the row cluster results from each submatrix in (2a) to construct
$\widehat{A}$

  5. (3b)

Combine the col. cluster results from each submatrix in (2b) to construct
$\widehat{B}$

Algorithm 5 Mean and Covariance Based Nested Clustering with COD Figure 5: The
$40\times 100$ matrix $M$. The first clustering layer ($\mathcal{G}^{(r),[1]}$
and $\mathcal{G}^{(c),[1]}$) is delineated in red while the second layer
($\mathcal{G}^{(r),[2]}$ and $\mathcal{G}^{(c),[2]}$) is delineated in green.
The dimensions of the ensuing grid like structure are denoted as numbers
outside the matrix. The numbers in red inside the matrix denote the mean
elements $M_{ij}$ which generate the data (and are the same throughout the
same “red” first layer nested cluster). For our simulation setting, the first
layer (red) is divided according to the mean value, while the second layer
(green) is divided according to the covariance structure of $Z$.

####  G.2.3 Simulation Study Comparing Algorithm 3 and Algorithm 5

In order to illustrate the feasibility of this nested algorithm, we include
simulation results with the aforementioned nested structure. For $n=5,\hskip
2.84544pt10,\hskip 2.84544pt15,\hskip 2.84544pt20,\hskip 2.84544pt25,\hskip
2.84544pt30,\hskip 2.84544pt35,\hskip 2.84544pt40,\hskip 2.84544pt45$, we
generate $n$ i.i.d. copies of a $40\times 100$ matrix
($p=40,\leavevmode\nobreak\ q=100$) with two nested row clusterings
$\\{\mathcal{G}^{(r),[1]},\mathcal{G}^{(r),[2]}\\}$ and two nested column
clusterings $\\{\mathcal{G}^{(c),[1]},\mathcal{G}^{(c),[2]}\\}$ where the row
clusterings are given by

$\displaystyle\mathcal{G}^{(r),[1]}$ $\displaystyle\leavevmode\nobreak\
=\leavevmode\nobreak\ \\{G^{(r),[1]}_{1}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \hskip 91.04872pt,\leavevmode\nobreak\
G^{(r),[1]}_{2}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \hskip 91.04872pt\\}$
$\displaystyle\mathcal{G}^{(r),[2]}$ $\displaystyle\leavevmode\nobreak\
=\leavevmode\nobreak\ \\{G^{(r),[2]}_{1},\leavevmode\nobreak\
G^{(r),[2]}_{2},\leavevmode\nobreak\ G^{(r),[2]}_{3},\leavevmode\nobreak\
G^{(r),[2]}_{4},\leavevmode\nobreak\ G^{(r),[2]}_{5},\leavevmode\nobreak\
G^{(r),[2]}_{6},\leavevmode\nobreak\ G^{(r),[2]}_{7},\leavevmode\nobreak\
G^{(r),[2]}_{8}\\},\leavevmode\nobreak\ \leavevmode\nobreak\
\textrm{where}\leavevmode\nobreak\ $ $\displaystyle G$
${}^{(r),[1]}_{1}=\\{1,2,\leavevmode\nobreak\ ...\leavevmode\nobreak\
,20\\},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \hskip
77.67589ptG^{(r),[1]}_{2}=\\{21,22,\leavevmode\nobreak\
...\leavevmode\nobreak\ ,40\\}$ $\displaystyle G$
${}^{(r),[2]}_{1}=\\{1,2,3\\},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \hskip
2.84544ptG^{(r),[2]}_{2}=\\{4,5,6,7\\},\leavevmode\nobreak\
\leavevmode\nobreak\ G^{(r),[2]}_{5}=\\{21,22,23\\},\leavevmode\nobreak\
G^{(r),[2]}_{6}=\\{24,25,26,27\\}$ $\displaystyle G$
${}^{(r),[2]}_{3}=\\{8,...,12\\},\leavevmode\nobreak\
G^{(r),[2]}_{4}=\\{13,...,20\\},\leavevmode\nobreak\ \hskip
0.85355ptG^{(r),[2]}_{7}=\\{28,...,32\\},\leavevmode\nobreak\
G^{(r),[2]}_{8}=\\{33,...,40\\}.$

The column clustering can be defined in a similar fashion, but for brevity it
is omitted here. The true clustering can be visualized with Figure 5.

Figure 6: The performance of the nested clustering algorithm (using both mean
and covariance information) and the 2-STEP COD method (using covariance
information only). The simulation setting is $p=40$, $q=100$,
$K_{1}=8,\leavevmode\nobreak\ K_{2}=20$,
$n\in\\{5,10,15,20,25,30,35,40,45\\}$. The ARI values for both methods are
averaged over 30 repetitions.

We use the ARI (Adjusted Rand Index, Vinh et al. (2009)) to compare the
performance of the nested clustering method with our original 2-STEP COD
method. The results are presented in Figure 6. Overall, considering the
feature matrices are very large ($40\times 100$), both methods work very well,
as we get near perfect recovery once $n\geq 45$. Taking a closer look, for
column clustering, the nested method is uniformly better than the original
2-STEP COD method over all $n$ values. For row clustering, for very small $n$
values of $5,10,15$ and $20$, the nested method performs better, but for
slightly larger $n$ values of $25,30,35,40$ and $45$, the original 2-STEP COD
method performs better. We observe that the closer the matrix is to a square,
the easier it is for 2-STEP COD to cluster the rows and columns. If the matrix
is highly unbalanced ($p>>q$ or $p<<q$), then due to the iterative nature of
2-STEP COD (using rows to cluster columns, then using columns to cluster the
rows, etc.), clustering one of the two dimensions becomes extremely difficult.
This may lead to the final result being inaccurate. In our simulation
settings, for row clustering, 2-STEP COD clusters over the entire $40\times
100$ matrix while the nested clustering method has to cluster over the two
$20\times 100$ submatrices separately (if the mean-based clustering can
perfectly recover $\mathcal{G}^{(r),[1]}$). This explains why the nested
algorithm may perform worse than the original 2-STEP COD for row clustering.

We also implement the competing methods DEEM, TGMM and TEMM (Mai et al.
(2022), Deng and Zhang (2022)), which are also included in the main simulation
results, but none are able to be implemented successfully in this nested
clustering simulation study. As explained in the main paper, their algorithms
depend on different modeling assumptions and may fail to converge in our
simulation setting.

All in all, we conclude that if the user wishes to incorporate mean
information and also look at the finer clustering for the rows and columns
that cannot be discerned by the mean alone, our COD based method can be used
as the second step in a nested clustering type process (Algorithm 5) to obtain
reliable results.

###  G.3 Higher Order Tensor Models

####  G.3.1 The Tensor Cluster Model and Algorithm 6

We present in the following the analog of our matrix clustering setting in a
setting that is one order higher. Consider a three-way tensor
$X\in\mathbb{R}^{J\times P\times Q}$. Assume that $X$ can be decomposed as

$X\leavevmode\nobreak\ =\leavevmode\nobreak\
Z\times_{1}A\times_{2}B\times_{3}C\leavevmode\nobreak\ +\leavevmode\nobreak\
\Gamma,$ (G.5)

where $Z\in\mathbb{R}^{K_{1}\times K_{2}\times K_{3}}$ is a three-way latent
tensor (i.e., each entry of $Z$ is a random variable) with $\mathbb{E}(Z)=0$,
$A\in\mathbb{R}^{J\times K_{1}}$, $B\in\mathbb{R}^{P\times K_{2}}$ and
$C\in\mathbb{R}^{Q\times K_{3}}$ are the unknown membership matrices for the
rows, columns and tubes of $X$, respectively, and
$\Gamma\in\mathbb{R}^{J\times P\times Q}$ represents the mean 0 random noise
tensor. Here, we use the $n$-mode product notation from the tensor literature,
e.g., $Z\times_{1}A$ is a $J\times K_{2}\times K_{3}$ tensor with
$(Z\times_{1}A)_{jst}=\sum_{k=1}^{K_{1}}Z_{kst}A_{jk}$. Similar to the matrix
clustering setting, we assume the entries of $\Gamma$ are mutually independent
and are also independent of $Z$. Each entry of the membership matrix $A$ takes
values in $\\{0,1\\}$, such that $A_{ak}=1$ if row $a$ belongs to row cluster
$k$ and $A_{ak}=0$ otherwise. The membership matrices $B$ and $C$ are defined
similarly. Assuming that $n$ i.i.d. samples $X^{(1)},...,X^{(n)}$ of a random
tensor $X\in\mathbb{R}^{J\times P\times Q}$ are observed, our goal is to
recover the membership matrices $A,B$ and $C$.

Next, we define the so-called matricization (or equivalently, “unfolding”) of
the tensor $X$. Let $X_{(1)}\in\mathbb{R}^{J\times PQ}$ denote the mode-1
matricization of $X$, which arranges the mode-1 fibers of $X$ as columns into
a matrix. Similarly, we can define $X_{(2)}\in\mathbb{R}^{P\times JQ}$ and
$X_{(3)}\in\mathbb{R}^{Q\times JP}$ as the mode-2 and mode-3 matricizations of
$X$. Similar to the matrix clustering setting, we consider the weighted
covariance matrix $\Sigma_{1,W_{1}}=\mathbb{E}(X_{(1)}W_{1}X_{(1)}^{T})$,
where $W_{1}\in\mathbb{R}^{PQ\times PQ}$ is some positive semi-definite weight
matrix to be chosen. It is easily seen that $\Sigma_{1,W_{1}}$ exhibits the
clustering structure for the rows, which can be identified via our COD
algorithm. In particular, let
$\widehat{\Sigma}_{j,W_{j}}=\frac{1}{n}\sum_{i=1}^{n}X_{(j)}^{(i)}W_{j}X_{(j)}^{(i)T}$
for $j=1,2,3$. The tensor clustering algorithm with the identity weight is
shown in Algorithm 6. For simplicity, we only consider the identity weight.
The algorithm may further be improved by using the optimal weight as in our
one-step Algorithm 2 or our two-step Algorithm 3, but we leave the detailed
analysis for future study.

  * for $j=1,2,3$

    * -

Compute $\widehat{\Sigma}_{j,W_{j}}$, where $W_{1}=\frac{1}{PQ}\cdot I_{PQ}$,
$W_{2}=\frac{1}{JQ}\cdot I_{JQ}$ and $W_{3}=\frac{1}{PJ}\cdot I_{PJ}$.

    * -

Apply Algorithm 1 to $\widehat{\Sigma}_{j,W_{j}}$ to construct
$\widehat{A},\widehat{B}$, and $\widehat{C}$, respectively.

Algorithm 6 Tensor Clustering Algorithm with the Identity Weight

####  G.3.2 Simulation Study

To illustrate the feasibility of our proposed method, we also conduct a
simulation study with Algorithm 6 for which the results are presented in Table
2. We have $n$ i.i.d. copies of a $15\times 10\times 10$ tensor
($n=10,20,30,40,50$) that have Toeplitz covariance matrices in each direction
with decay rates of -0.4, 0.3 and -0.2, respectively. The noise is generated
in the homogeneous setting with the noise variance being $15$ to match our
main simulation settings. It is seen that clustering is more accurate as the
sample size grows, which provides empirical evidence for applying Algorithm 6
to higher-order tensors.

|  |  Row (ARI) |  Column (ARI) |  Tube (ARI)
---|---|---|---|---
$n$ | 10 | 0.0863 | 0 | 0.4449
20 | 0.6637 | 0.3095 | 0.8904
| 30 | 0.7594 | 0.7936 | 0.9334
| 40 | 0.7864 | 0.8680 | 0.9033
| 50 | 0.8107 | 0.8753 | 0.9889
Table 2: The ARI values obtained from the Tensor Clustering Algorithm with the
Identity Weight (Algorithm 6) under the $n=10,20,30,40,50$ setting.

While our preliminary simulation results are indeed promising, proving minimax
optimality would be challenging in the tensor clustering. Also, as mentioned
above, since using more intricate weight matrices instead of the identity
matrix may improve the clustering performance, a rigorous analysis of the
algorithm would also be interesting. All of these directions are promising
fields of future study.

##  Appendix H Proofs

###  H.1 Clustering consistency

First, we state the Hanson-Wright inequality, which is instrumental in the
proof.

######  Lemma H.1.

(Hanson-Wright) There exist positive constants $c,c^{\prime}$ such that for
all $n\times n$ matrices $H$, if $\xi=(\xi_{1},...,\xi_{n})^{T}$ is a vector
of independent mean 0 sub-Gaussian random variables with
$||\xi_{i}||_{\psi_{2}}\leq L$ for some $L>0$, then for all t, the following
holds:

$\displaystyle\mathbb{P}\bigg{[}\big{|}\xi^{T}H\xi-\mathbb{E}(\xi^{T}H\xi)\big{|}>cL^{2}(||H||_{F}\sqrt{t}+c^{\prime}||H||_{\text{op}}t)\bigg{]}\leq
e^{-t}.$

####  H.1.1 Proof of Theorem 3.1

###### Proof.

Denote $\widehat{L}\widehat{L}^{T}=\widehat{W}$ and
$\bar{X}^{(i)}=X^{(i)}\widehat{L}$. So,
$\widehat{\Sigma}_{p,\widehat{W}}=\frac{1}{|D_{2}|}\sum_{i\in
D_{2}}\bar{X}^{(i)}\bar{X}^{(i)T}$. For simplicity, we assume $|D_{2}|=n/2$,
and the summation is from $1$ to $n/2$. Furthermore, denote
$H_{k,m}^{+}=\operatorname{{\rm
Var}}(\bar{X}_{k\cdot}+\bar{X}_{m\cdot}|\widehat{W})$ and
$H_{k,m}^{-}=\operatorname{{\rm
Var}}(\bar{X}_{k\cdot}-\bar{X}_{m\cdot}|\widehat{W})$. Denote
$\widetilde{X}_{k,m}^{+}=(H_{k,m}^{+})^{-1/2}(\bar{X}_{k\cdot}^{T}+\bar{X}_{m\cdot}^{T})$.
Since $\widetilde{X}_{k,m}^{+}$ has mean 0 variance 1 given $\widehat{W}$ and
$X$ is Gaussian, $\widetilde{X}_{k,m}^{+}$ is a vector of independent Gaussian
random variables given $\widehat{W}$.

Denote

$\displaystyle\widetilde{\mathbb{X}}_{k,m}^{+}=\begin{bmatrix}(H_{k,m}^{+})^{-1/2}&&\\\
&...&\\\ &&(H_{k,m}^{+})^{-1/2}\\\
\end{bmatrix}\begin{bmatrix}(\bar{X}_{k\cdot}^{(1)}+\bar{X}_{m\cdot}^{(1)})^{T}\\\
:\\\ (\bar{X}_{k\cdot}^{(n/2)}+\bar{X}_{m\cdot}^{(n/2)})^{T}\\\
\end{bmatrix}=\begin{bmatrix}\widetilde{X}_{k,m}^{+(1)}\\\ :\\\
\widetilde{X}_{k,m}^{+(n/2)}\\\ \end{bmatrix}$
$\displaystyle\widetilde{\mathbb{X}}_{k,m}^{-}=\begin{bmatrix}(H_{k,m}^{-})^{-1/2}&&\\\
&...&\\\ &&(H_{k,m}^{-})^{-1/2}\\\
\end{bmatrix}\begin{bmatrix}(\bar{X}_{k\cdot}^{(1)}-\bar{X}_{m\cdot}^{(1)})^{T}\\\
:\\\ (\bar{X}_{k\cdot}^{(n/2)}-\bar{X}_{m\cdot}^{(n/2)})^{T}\\\
\end{bmatrix}=\begin{bmatrix}\widetilde{X}_{k,m}^{-(1)}\\\ :\\\
\widetilde{X}_{k,m}^{-(n/2)}\\\ \end{bmatrix}$ and
$\displaystyle\mathbb{H}_{k,m}^{+}=\begin{bmatrix}(H_{k,m}^{+})^{-1/2}&&\\\
&...&\\\ &&(H_{k,m}^{+})^{-1/2}\\\ \end{bmatrix},\hskip
2.84544pt\mathbb{H}_{k,m}^{-}=\begin{bmatrix}(H_{k,m}^{-})^{-1/2}&&\\\
&...&\\\ &&(H_{k,m}^{-})^{-1/2}\\\ \end{bmatrix},$

where $\mathbb{H}_{k,m}^{+}$ and $\mathbb{H}_{k,m}^{-}$ each have $n/2$ blocks
on the diagonal. Note that the quadratic forms can be used to get the sum of
the products. The following calculation holds for one observation:

$\displaystyle(\widetilde{X}_{k,m}^{+})^{T}H_{k,m}^{+}\widetilde{X}_{k,m}^{+}-(\widetilde{X}_{k,m}^{-})^{T}H_{k,m}^{-}\widetilde{X}_{k,m}$
$\displaystyle=\sum_{j=1}^{q}(\bar{X}_{kj}^{2}+2\bar{X}_{kj}\bar{X}_{mj}+\bar{X}_{mj}^{2})$
$\displaystyle\hskip
85.35826pt-\sum_{j=1}^{q}(\bar{X}_{kj}^{2}-2\bar{X}_{kj}\bar{X}_{mj}+\bar{X}_{mj}^{2})$
$\displaystyle=4\sum_{j=1}^{q}\bar{X}_{kj}\bar{X}_{mj}.$

We can extend this to $n/2$ observations by using block matrix notation from
above.

$\displaystyle(\widetilde{\mathbb{X}}_{k,m}^{+})^{T}\mathbb{H}_{k,m}^{+}\widetilde{\mathbb{X}}_{k,m}^{+}-(\widetilde{\mathbb{X}}_{k,m}^{-})^{T}\mathbb{H}_{k,m}^{-}\widetilde{\mathbb{X}}_{k,m}^{-}$
$\displaystyle=\sum_{i=1}^{n/2}\sum_{j=1}^{q}(\bar{X}_{kj}^{(i)2}+2\bar{X}_{kj}^{(i)}\bar{X}_{mj}^{(i)}+\bar{X}_{mj}^{(i)2})-\sum_{i=1}^{n/2}\sum_{j=1}^{q}(\bar{X}_{kj}^{(i)2}-2\bar{X}_{kj}^{(i)}\bar{X}_{mj}^{(i)}+\bar{X}_{mj}^{(i)2})$
$\displaystyle=4\sum_{i=1}^{n/2}\sum_{j=1}^{q}\bar{X}_{kj}^{(i)}\bar{X}_{mj}^{(i)}.$

Thus, we have

$\displaystyle||\widehat{\Sigma}_{p,\widehat{W}}-\Sigma_{p,\widehat{W}}||_{\max}$
$\displaystyle=\underset{k,m}{\max}\bigg{|}\frac{2}{n}\sum_{i=1}^{n/2}\sum_{j=1}^{q}\bar{X}_{kj}^{(i)}\bar{X}_{mj}^{(i)}-\sum_{j=1}^{q}\mathbb{E}(\bar{X}_{kj}\bar{X}_{mj}|\widehat{W})\bigg{|}$
$\displaystyle=\underset{k,m}{\text{max}}\bigg{|}\frac{1}{2n}\bigg{(}(\widetilde{\mathbb{X}}_{k,m}^{+})^{T}\mathbb{H}_{k,m}^{+}\widetilde{\mathbb{X}}_{k,m}^{+}-\mathbb{E}((\widetilde{\mathbb{X}}_{k,m}^{+})^{T}\mathbb{H}_{k,m}^{+}\widetilde{\mathbb{X}}_{k,m}^{+}|\widehat{W})\bigg{)}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-\frac{1}{2n}\bigg{(}(\widetilde{\mathbb{X}}_{k,m}^{-})^{T}\mathbb{H}_{k,m}^{-}\widetilde{\mathbb{X}}_{k,m}^{-}-\mathbb{E}((\widetilde{\mathbb{X}}_{k,m}^{-})^{T}\mathbb{H}_{k,m}^{-}\widetilde{\mathbb{X}}_{k,m}^{-}|\widehat{W})\bigg{)}\bigg{|}.$

By Lemma H.1, we know that there exist positive constants
$c,c^{\prime},c^{\prime\prime},c^{\prime\prime\prime}$ such that: conditioned
on $\widehat{W}$ with probability greater than $1-e^{-t}$,

$\displaystyle\big{|}(\widetilde{\mathbb{X}}_{k,m}^{+})^{T}\mathbb{H}_{k,m}^{+}\widetilde{\mathbb{X}}_{k,m}^{+}-\mathbb{E}((\widetilde{\mathbb{X}}_{k,m}^{+})^{T}\mathbb{H}_{k,m}^{+}\widetilde{\mathbb{X}}_{k,m}^{+}|\widehat{W})\big{|}$
$\displaystyle\leq
c\bigg{(}||\mathbb{H}_{k,m}^{+}||_{F}\sqrt{t}+c^{\prime}||\mathbb{H}_{k,m}^{+}||_{\text{op}}\cdot
t\bigg{)}.$

Similarly, conditioned on $\widehat{W}$ with probability greater than
$1-e^{-t}$,

$\big{|}(\widetilde{\mathbb{X}}_{k,m}^{-})^{T}\mathbb{H}_{k,m}^{-}\widetilde{\mathbb{X}}_{k,m}^{-}-\mathbb{E}((\widetilde{\mathbb{X}}_{k,m}^{-})^{T}\mathbb{H}_{k,m}^{-}\widetilde{\mathbb{X}}_{k,m}^{-}|\widehat{W})\big{|}\leq
c^{\prime\prime}\bigg{(}||\mathbb{H}_{k,m}^{-}||_{F}\sqrt{t}+c^{\prime\prime\prime}||\mathbb{H}_{k,m}^{-}||_{\text{op}}\cdot
t\bigg{)}.$

By the definition of $\mathbb{H}_{k,m}^{+}$, we know
$||\mathbb{H}_{k,m}^{+}||_{F}^{2}=n||H_{k,m}^{+}||_{F}^{2}/2$ and
$||\mathbb{H}_{k,m}^{+}||_{\text{op}}=||H_{k,m}^{+}||_{\text{op}}$. Thus, with
a union bound over $k,m$, we get the following concentration inequality:

$\displaystyle||\widehat{\Sigma}_{p,\widehat{W}}-\Sigma_{p,\widehat{W}}||_{\max}$
$\displaystyle\lesssim\frac{1}{n}\cdot\max_{k,m}\bigg{(}\sqrt{||H_{k,m}^{+}||_{F}^{2}\cdot
nt}+||H_{k,m}^{+}||_{\text{op}}\cdot t\bigg{)}$ $\displaystyle\hskip
56.9055pt+\frac{1}{n}\cdot\max_{k,m}\bigg{(}\sqrt{||H_{k,m}^{-}||_{F}^{2}\cdot
nt}+||H_{k,m}^{-}||_{\text{op}}\cdot t\bigg{)}$ (H.1)

conditioned on $\widehat{W}$ with probability greater than
$1-\frac{p(p-1)}{2}e^{-t}$, and setting $t=3\log p$, we know the inequality
(H.1) holds with conditional probability greater than $1-\frac{1}{2p}$. By
taking another expectation with $\widehat{W}$, the inequality (H.1) holds with
probability greater than $1-\frac{1}{2p}$ unconditionally.

Because $||\cdot||_{F}\geq||\cdot||_{op}$ holds in general, the Frobenius term
dominates in the regime of $\frac{\log p}{n}=o(1)$, so all that is left to do
is to bound $||H_{k,m}^{+}||_{F}^{2}$. An identical argument can be made for
$||H_{k,m}^{-}||_{F}^{2}$, and we know the following upper bound holds for
both terms,

$\displaystyle||H_{k,m}^{+}||_{F}$ $\displaystyle\leq||\operatorname{{\rm