Multiway clustering via tensor block models

Miaoyan Wang and Yuchen Zeng111Miaoyan Wang is Assistant Professor, Department
of Statistics, University of Wisconsin-Madison, Madison, WI 53706, E-mail:
<EMAIL_ADDRESS>Yuchen Zheng is a BS/MS student, Department of
Statistics, University of Wisconsin-Madison, Madison, WI 53706, E-mail:
<EMAIL_ADDRESS>
Department of Statistics, University of Wisconsin-Madison

###### Abstract

We consider the problem of identifying multiway block structure from a large
noisy tensor. Such problems arise frequently in applications such as genomics,
recommendation system, topic modeling, and sensor network localization. We
propose a tensor block model, develop a unified least-square estimation, and
obtain the theoretical accuracy guarantees for multiway clustering. The
statistical convergence of the estimator is established, and we show that the
associated clustering procedure achieves partition consistency. A sparse
regularization is further developed for identifying important blocks with
elevated means. The proposal handles a broad range of data types, including
binary, continuous, and hybrid observations. Through simulation and
application to two real datasets, we demonstrate the outperformance of our
approach over previous methods.

Keywords: Tensor block model, Clustering, Least-square estimation, Dimension
reduction.

##  1 Introduction

Higher-order tensors have recently attracted increased attention in data-
intensive fields such as neuroscience [1], social networks [2], computer
vision [3], and genomics [4, 5]. In many applications, the data tensors are
often expected to have underlying block structure. One example is multi-tissue
expression data [4], in which genome-wide expression profiles are collected
from different tissues in a number of individuals. There may be groups of
genes similarly expressed in subsets of tissues and individuals;
mathematically, this implies an underlying three-way block structure in the
data tensor. In a different context, block structure may emerge in a binary-
valued tensor. Examples include multilayer network data [2], with the nodes
representing the individuals and the layers representing the multiple types of
relations. Here a planted block represents a community of individuals that are
highly connected within a class of relationships.

Figure 1: Examples of tensor block model (TBM). (a) Our TBM method is used for
multiway clustering and for revealing the underlying checkerbox structure in a
noisy tensor. (b) The sparse TBM method is used for detecting sub-tensors of
elevated means.

This paper presents a new method and the associated theory for tensors with
block structure. We develop a unified least-square estimation procedure for
identifying multiway block structure. The proposal applies to a broad range of
data types, including binary, continuous, and hybrid observations. We
establish a high-probability error bound for the resulting estimator, and show
that the procedure enjoys consistency guarantees on the block structure
recovery as the dimension of the data tensor grows. Furthermore, we develop a
sparse extension of the tensor block model for block selections. Figure 1
shows two immediate examples of our method. When the data tensor possesses a
checkerbox pattern modulo some unknown reordering of entries, our method
amounts to multiway clustering that simultaneously clusters each mode of the
tensor (Figure 1a). When the data tensor has no full checkerbox structure but
contains a small numbers of sub-tensors of elevated means, we develop a sparse
version of our method to detect these sub-tensors of interest (Figure 1b).

Related work. Our work is closely related to, but also clearly distinctive
from, the low-rank tensor decomposition. A number of methods have been
developed for low-rank tensor estimation, including CANDECOMP/PARAFAC (CP)
decomposition [6] and Tucker decomposition [7]. The CP model decomposes a
tensor into a sum of rank-1 tensors, whereas Tucker model decomposes a tensor
into a core tensor multiplied by orthogonal matrices in each mode. In this
paper we investigate an alternative block structure assumption, which has yet
to be studied for higher-order tensors. Note that a block structure
automatically implies low-rankness. However, as we will show in Section 4, a
direct application of low rank estimation to the current setting will result
in an inferior estimator. Therefore, a full exploitation of the block
structure is necessary; this is the focus of the current paper.

Our work is also connected to biclustering [8] and its higher-order extensions
[9, 10]. Existing multiway clustering methods [9, 10, 5, 11] typically take a
two-step procedure, by first estimating a low-dimension representation of the
data tensor and then applying clustering algorithms to the tensor factors. In
contrast, our tensor block model takes a single shot to perform estimation and
clustering simultaneously. This approach achieves a higher accuracy and an
improved interpretability. Moreover, earlier solutions to multiway clustering
[12, 9] focus on the algorithm effectiveness, leaving the statistical
optimality of the estimators unaddressed. Very recently, Chi et al [13]
provides an attempt to study the statistical properties of the tensor block
model. We will show that our estimator obtains a faster convergence rate than
theirs, and the power is further boosted with a sparse regularity.

##  2 Preliminaries

We begin by reviewing a few basic factors about tensors [14]. We use
$\mathcal{Y}=\llbracket
y_{i_{1},\ldots,i_{K}}\rrbracket\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$
to denote an order-$K$ $(d_{1},\ldots,d_{K})$-dimensional tensor. The
multilinear multiplication of a tensor
$\mathcal{Y}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$ by matrices
$\bm{M}_{k}=\llbracket
m_{i_{k},j_{k}}^{(k)}\rrbracket\in\mathbb{R}^{s_{k}\times d_{k}}$ is defined
as

$\mathcal{Y}\times_{1}\bm{M}_{1}\ldots\times_{K}\bm{M}_{K}=\llbracket\sum_{i_{1},\ldots,i_{K}}y_{i_{1},\ldots,i_{K}}m_{i_{1},j_{1}}^{(1)}\ldots
m_{i_{K},j_{K}}^{(K)}\rrbracket,$

which results in an order-$K$ tensor $(s_{1},\ldots,s_{K})$-dimensional
tensor. For any two tensors $\mathcal{Y}=\llbracket
y_{i_{1},\ldots,i_{K}}\rrbracket$, $\mathcal{Y}^{\prime}=\llbracket
y^{\prime}_{i_{1},\ldots,i_{K}}\rrbracket$ of identical order and dimensions,
their inner product is defined as
$\langle\mathcal{Y},\mathcal{Y}^{\prime}\rangle=\sum_{i_{1},\ldots,i_{K}}y_{i_{1},\ldots,i_{K}}y^{\prime}_{i_{1},\ldots,i_{K}}$.
The Frobenius norm of tensor $\mathcal{Y}$ is defined as
$\lVert\mathcal{Y}\rVert_{F}=\langle\mathcal{Y},\mathcal{Y}\rangle^{1/2}$; it
is the Euclidean norm of $\mathcal{Y}$ regarded as an
$\prod_{k}d_{k}$-dimensional vector. The maximum norm of tensor $\mathcal{Y}$
is defined as
$\lVert\mathcal{Y}\rVert_{\max}=\max_{i_{1},\ldots,i_{K}}|y_{i_{1},\ldots,i_{K}}|$.
An order-($K$-$1$) slice of $\mathcal{Y}$ is a sub-tensor of $\mathcal{Y}$
obtained by holding the index in one mode fixed while letting other indices
vary.

A clustering of $d$ objects is a partition of the index set
$[d]:=\\{1,2,\ldots,d\\}$ into $R$ disjoint non-empty subsets. We refer to the
number of clusters, $R$, as the clustering size. Equivalently, the clustering
(or partition) can be represented using the “membership matrix”. A membership
matrix $\bm{M}\in\mathbb{R}^{R\times d}$ is an incidence matrix whose
$(i,j)$-entry is 1 if and only if the element $j$ belongs to the cluster $i$,
and 0 otherwise. Throughout the paper, we will use the terms “clustering”,
“partition”, and “membership matrix” exchangeably. For a higher-order tensor,
the concept of index partition applies to each of the modes. A block is a sub-
tensor induced by the index partitions along each of the $K$ modes. We use the
term “cluster” to refer to the marginal partition on mode $k$, and reserve the
term “block” for the multiway partition of the tensor.

##  3 Tensor block model

Let $\mathcal{Y}=\llbracket
y_{i_{1},\ldots,i_{K}}\rrbracket\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$
denote an order-$K$, $(d_{1},\ldots,d_{K})$-dimensional data tensor. The main
assumption of tensor block model (TBM) is that the observed data tensor
$\mathcal{Y}$ is a noisy realization of an underlying tensor that exhibits a
checkerbox structure (see Figure 1a). Specifically, suppose that the $k$-th
mode of the tensor consists of $R_{k}$ clusters. If the tensor entry
$y_{i_{1},\ldots,i_{K}}$ belongs to the block determined by the $r_{k}$th
cluster in the mode $k$ for $r_{k}\in[R_{k}]$, then we assume that

$y_{i_{1},\ldots,i_{K}}=c_{r_{1},\ldots,r_{K}}+\varepsilon_{i_{1},\ldots,i_{K}},\quad\text{for
}(i_{1},\ldots,i_{K})\in[d_{1}]\times\cdots\times[d_{K}],$ (1)

where $c_{r_{1},\ldots,r_{K}}$ is the mean of the tensor block indexed by
$(r_{1},\ldots,r_{K})$, and $\varepsilon_{i_{1},\ldots,i_{K}}$’s are
independent, mean-zero noise terms to be specified later. Our goal is to (i)
find the clustering along each of the modes, and (ii) estimate the block means
$\\{c_{r_{1},\ldots,r_{K}}\\}$, such that a corresponding blockwise-constant
checkerbox structure emerges in the data tensor.

The tensor block model (1) falls into a general class of non-overlapping,
constant-mean clustering models [15], in that each tensor entry belongs to
exactly one block with a common mean. The TBM can be equivalently expressed as
a special tensor Tucker model,

$\mathcal{Y}=\mathcal{C}\times_{1}\bm{M}_{1}\times_{2}\cdots\times_{K}\bm{M}_{K}+\mathcal{E},$
(2)

where $\mathcal{C}=\llbracket
c_{r_{1},\ldots,r_{K}}\rrbracket\in\mathbb{R}^{R_{1}\times\cdots\times R_{K}}$
is a core tensor consisting of block means,
$\bm{M}_{k}\in\\{0,1\\}^{d_{k}\times R_{k}}$ is a membership matrix indicating
the block allocations along mode $k$ for $k\in[K]$, and
$\mathcal{E}=\llbracket\varepsilon_{i_{1},\ldots,i_{K}}\rrbracket$ is the
noise tensor. We view the TBM (2) as a super-sparse Tucker model, in the sense
that the each column of $\bm{M}_{k}$ consists of one copy of 1’s and massive
0’s.

We make a general assumption on the noise tensor $\mathcal{E}$. The noise
terms $\varepsilon_{i_{1},\ldots,i_{K}}$’s are assumed to be independent,
mean-zero $\sigma$-subgaussian, where $\sigma>0$ is the subgaussianity
parameter. More precisely,

$\mathbb{E}e^{\lambda\varepsilon_{i_{1},\ldots,i_{K}}}\leq
e^{\lambda^{2}\sigma^{2}/2},\quad\text{for all
}(i_{1},\ldots,i_{K})\in[d_{1}]\times\cdots\times[d_{K}]\ \text{and all}\
\lambda\in\mathbb{R}.$ (3)

Th assumption (3) incorporates common situations such as Gaussian noise,
Bernoulli noise, and noise with bounded support. In particular, we consider
two important examples of the TBM:

######  Example 1 (Gaussian tensor block model).

Let $\mathcal{Y}$ be a continuous-valued tensor. The Gaussian tensor block
model (GTBM)
$y_{i_{1},\ldots,i_{K}}\sim_{\text{i.i.d.}}N(c_{r_{1},\ldots,r_{K}},\sigma^{2})$
is a special case of model (1), with the subgaussianity parameter $\sigma$
equal to the error variance. The GTBM serves as the foundation for many tensor
clustering algorithms [12, 4, 13].

######  Example 2 (Stochastic tensor block model).

Let $\mathcal{Y}$ be a binary-valued tensor. The stochastic tensor block model
(STBM)
$y_{i_{1},\ldots,i_{K}}\sim_{\text{i.i.d.}}\text{Bernoulli}(c_{r_{1},\ldots,r_{K}})$
is a special case of model (1), with the subgaussianity parameter $\sigma$
equal to ${1\over 4}$. The STBM can be viewed as an extension, to higher-order
tensors, of the popular stochastic block model [16, 17] for matrix-based
network analysis. In the filed of community detection, multi-layer stochastic
model has also been developed for multi-relational network data analysis [18,
19].

More generally, our model also applied to hybrid error distributions, in which
different types of distribution are allowed for different portions of the
tensor. This scenario may happen, for example, when the data tensor
$\mathcal{Y}$ represents concatenated measurements from multiple data sources.

Before we discuss the estimation, we present the identifiability of the TBM.

######  Assumption 1 (Irreducible core).

The core tensor $\mathcal{C}$ is called irreducible if it cannot be written as
a block tensor with the number of mode-$k$ clusters smaller than $R_{k}$, for
any $k\in[K]$.

In the matrix case $(K=2)$, the irreducibility is equivalent to saying that
$\mathcal{C}$ has no two identical rows and no two identical columns. In the
higher-order case, the assumption requires that none of order-($K$-$1$) slices
of $\mathcal{C}$ are identical. Note that irreducibility is a weaker
assumption than full-rankness.

######  Proposition 1 (Identifiability).

Consider a Gaussian or Bernoulli TBM (1). Under Assumption 1, the factor
matrices $\bm{M}_{k}$’s are identifiable up to permutations of cluster labels.

The identifiability property for the TBM outperforms that for the classical
factor model [20, 21]. In the Tucker [22, 14] and many other factor analyses
[20, 21], the factors are identifiable only up to orthogonal rotations. Those
models recover only the (column) space spanned by $\bm{M}_{k}$, but not the
individual factors. In contrast, our model does not suffer from rotational
invariance, and as we show in Section 4, every individual factor is
consistently estimated in high dimensions. This brings a benefit to the
interpretation of factors in the tensor block model.

We propose a least-square approach for estimating the TBM. Let
$\Theta=\mathcal{C}\times_{1}\bm{M}_{1}\times_{2}\cdots\times_{K}\bm{M}_{K}$
denote the mean signal tensor with block structure. The mean tensor is assumed
to belong to the following parameter space

$\displaystyle\mathcal{P}_{R_{1},\ldots,R_{K}}=$ $\displaystyle\
\big{\\{}\Theta\in\mathbb{R}^{d_{1}\times\cdots\times
d_{K}}\colon\Theta=\mathcal{C}\times_{1}\bm{M}_{1}\times_{2}\cdots\times_{K}\bm{M}_{K},\text{with
some}$ $\displaystyle\quad\text{membership matrices $\bm{M}_{k}$'s and a core
tensor $\mathcal{C}\in\mathbb{R}^{R_{1}\times\cdots\times R_{K}}$}\big{\\}}.$

In the following theoretical analysis, we assume the clustering size
$\bm{R}=(R_{1},\ldots,R_{K})$ is known and simply write $\mathcal{P}$ for
short. The adaptation of unknown $\bm{R}$ will be addressed in Section 5.2.
The least-square estimator for the TBM (1) is

$\hat{\Theta}=\operatorname*{arg\,min}_{\Theta\in\mathcal{P}}\left\\{-2\langle\mathcal{Y},\Theta\rangle+\lVert\Theta\rVert_{F}^{2}\right\\}.$
(4)

The objective is equal (ignoring constants) to the sum of squares
$\lVert\mathcal{Y}-\Theta\rVert_{F}^{2}$ and hence the name of our estimator.

##  4 Statistical convergence

In this section, we establish the convergence rate of the least-squares
estimator (4) for two measurements. The first measurement is mean squared
error (MSE):

$\text{MSE}(\Theta_{\text{true}},\hat{\Theta})={1\over\prod_{k}d_{k}}\lVert\Theta_{\text{true}}-\hat{\Theta}\rVert_{F}^{2},$
(5)

where $\Theta_{\text{true}},\hat{\Theta}\in\mathcal{P}$ are the true and
estimated mean tensors, respectively. While the loss function corresponds to
the likelihood for the Gaussian tensor model, the same assertion does not hold
for other types of distribution such as stochastic tensor block model. We will
show that, with very high probability, a simple least-square estimator
achieves a fast convergence rate in a general class of block tensor models.

######  Theorem 1 (Convergence rate of MSE).

Let $\hat{\Theta}$ be the least-square estimator of $\Theta_{\text{true}}$
under model (1). There exists two constants $C_{1},C_{2}>0$ such that,

$\text{MSE}(\Theta_{\text{true}},\hat{\Theta})\leq{C_{1}\sigma^{2}\over\prod_{k}d_{k}}\left(\prod_{k}R_{k}+\sum_{k}d_{k}\log
R_{k}\right)$ (6)

holds with probability at least
$1-\exp(-C_{2}(\prod_{k}R_{k}+\sum_{k}d_{k}\log R_{k}))$ uniformly over
$\Theta_{\text{true}}\in\mathcal{P}$ and all error distribution satisfying
(3).

The convergence rate of MSE in (6) consists of two parts. The first part
$\prod_{k}R_{k}$ is the number of parameters in the core tensor $\mathcal{C}$,
while the second part $\sum_{k}d_{k}\log R_{k}$ reflects the the complexity
for estimating $\bm{M}_{k}$’s. It is the price that one has to pay for not
knowing the locations of the blocks.

We compare our bound with existing literature. The Tucker tensor decomposition
has a minimax convergence rate proportional to $\sum_{k}d_{k}R^{\prime}_{k}$
[22], where $R^{\prime}_{k}$ is the multilinear rank in the mode $k$. Applying
Tucker decomposition to the TBM yields $\sum_{k}d_{k}R_{k}$, because the
mode-$k$ rank is bounded by the number of mode-$k$ clusters. Now, as both the
dimension $d_{\min}=\min_{k}d_{k}$ and clustering size
$R_{\min}=\min_{k}R_{k}$ tend to infinity, we have
$\prod_{k}R_{k}+\sum_{k}d_{k}\log R_{k}\ll\sum_{k}d_{k}R_{k}$. Therefore, by
fully exploiting the block structure, we obtain a better convergence rate than
previously possible.

Recently, [13] proposed a convex relaxation for estimating the TBM. In the
special case when the tensor dimensions are equal at every mode
$d_{1}=\ldots=d_{K}=d$, their estimator has a convergence rate of order
$\mathcal{O}(d^{-1})$ for all $K\geq 2$. As we see from (6), our estimate
obtains a much better convergence rate $\mathcal{O}(d^{-(K-1)})$, which is
especially favorable as the order increases.

The bound (6) generalizes the previous results on structured matrix estimation
in network analysis [23, 16]. Earlier work [16] suggests the following
heuristics on the sample complexity for the matrix case:

${\text{(number of parameters)}+\log\text{(complexity of
models)}\over\text{number of samples}}.$ (7)

Our result supports this important principle for general $K\geq 2$. Note that,
in the TBM, the sample size is the total number of entries $\prod_{k}d_{k}$,
the number of parameters is $\prod_{k}R_{k}$, and the combinatoric complexity
for estimating block structure is of order $\prod_{k}R_{k}^{d_{k}}$.

Next we study the consistency of partition. To define the misclassification
rate (MCR), we need to introduce some additional notation. Let
$\bm{M}_{k}=\llbracket
m^{(k)}_{i,r}\rrbracket,\hat{\bm{M}}_{k}=\llbracket\hat{m}^{(k)}_{i,r^{\prime}}\rrbracket$
be two mode-$k$ membership matrices, and $\bm{D}^{(k)}=\llbracket
D_{r,r^{\prime}}^{(k)}\rrbracket$ be the mode-$k$ confusion matrix with
element
$D_{r,r^{\prime}}^{(k)}=\frac{1}{d_{k}}\sum_{i=1}^{d_{k}}\mathds{1}\\{m_{i,r}^{(k)}=\hat{m}_{i,r^{\prime}}^{(k)}=1\\}$,
where $r,r^{\prime}\in[R_{k}]$. Note that the row/column sum of $\bm{D}^{(k)}$
represents the nodes proportion in each cluster defined by $\bm{M}_{k}$ or
$\hat{\bm{M}}_{k}$. We restrict ourselves to non-degenerating clusterings;
that is, the row/column sums of $\bm{D}^{(k)}$ are lowered bounded by
$\tau>0$. With a little abuse of notation, we still use
$\mathcal{P}=\mathcal{P}(\tau)$ to denote the parameter space with the non-
degenerating assumption. The least-square estimator (4) should also be
interpreted with this constraint imposed.

We define the mode-$k$ misclassification rate (MCR) as

$\text{MCR}(\bm{M}_{k},\hat{\bm{M}}_{k})=\max_{r\in[R_{k}],a\neq
a^{\prime}\in[R_{k}]}\min\left\\{D^{(k)}_{a,r},\
D^{(k)}_{a^{\prime},r}\right\\}.$ (8)

In other words, MCR is the element-wise maximum of the confusion matrix after
removing the largest entry from each column. Under the non-degenerating
assumption, $\text{MCR}=0$ if and only if the confusion matrix $\bm{D}^{(k)}$
is a permutation of a diagonal matrix; that is, the estimated partition
matches with the true partition, up to permutations of cluster labels.

######  Theorem 2 (Convergence rate of MCR).

Consider a tensor block model (2) with sub-Gaussian parameter $\sigma$. Define
the minimal gap between the blocks $\delta_{\min}=\min_{k}\delta^{(k)}$, where
$\delta^{(k)}=\min_{r_{k}\neq r_{k}^{\prime}}$
$\max_{r_{1},\ldots,r_{k-1},r_{k+1},\ldots,r_{K}}(c_{r_{1},\ldots,r_{k},\ldots,r_{K}}-c_{r_{1},\ldots,r_{k}^{\prime},\ldots,r_{K}})^{2}$.
Let $\bm{M}_{k,\text{true}}$ be the true mode-$k$ membership,
$\hat{\bm{M}}_{k}$ be the estimator from (4). Then, for any
$\varepsilon\in[0,1]$,

$\mathbb{P}(\text{MCR}(\hat{\bm{M}}_{k},\bm{M}_{k,\text{true}})\geq\varepsilon)\leq
2^{1+\sum_{k}d_{k}}\exp\left(-\frac{C\varepsilon^{2}\delta_{\min}^{2}\tau^{3K-2}\prod_{k=1}^{K}d_{k}}{\sigma^{2}\lVert\mathcal{C}\rVert_{\max}^{2}}\right),$
(9)

where $C>0$ is a positive constant, and $\tau>0$ the lower bound of cluster
proportions.

The above theorem shows that our estimator consistently recovers the block
structure as the dimension of the data tensor grows. The block-mean gap
$\delta_{\min}$ serves the role of the eigen-separation as in the classical
tensor Tucker decomposition [22]. Table 1 summarizes the comparison of various
tensor methods in the special case when $d_{1}=\cdots=d_{K}=d$ and
$R_{1}=\cdots=R_{K}=R$.

Method | Recovery error (MSE) | Clustering error (MCR) | Block detection (see Section 6)
---|---|---|---
Tucker [22] | $dR$ | - | No
CoCo [13] | $d^{K-1}$ | - | No
TBM (this paper) | $d\log R$ | ${\sigma\lVert\mathcal{C}\rVert_{\max}\sqrt{\log d}\over\delta_{\min}\tau^{(3K-2)/2}}d^{-(K-1)/2}$ | Yes

Table 1: Comparison of various tensor decomposition methods.

##  5 Numerical implementation

###  5.1 Alternating optimization

We introduce an alternating optimization for solving (4). Estimating $\Theta$
consists of finding both the core tensor $\mathcal{C}$ and the membership
matrices $\bm{M}_{k}$’s. The optimization (4) can be written as

$\displaystyle(\hat{\mathcal{C}},\\{\hat{\bm{M}}_{k}\\})$
$\displaystyle=\operatorname*{arg\,min}_{\mathcal{C}\in\mathbb{R}^{R_{1}\times\cdots\times
R_{K}},\text{ membership matrices
$\bm{M}_{k}$'s}}f(\mathcal{C},\\{\bm{M}_{k}\\}),$ where $\displaystyle\quad
f(\mathcal{C},\\{\bm{M}_{k}\\})=\lVert\mathcal{Y}-\mathcal{C}\times_{1}\bm{M}_{1}\times_{2}\ldots\times_{K}\bm{M}_{K}\rVert_{F}^{2}.$

The decision variables consist of $K+1$ blocks of variables, one for the core
tensor $\mathcal{C}$ and $K$ for the membership matrices $\bm{M}_{k}$’s. We
notice that, if any $K$ out of the $K+1$ blocks of variables are known, then
the last block of variables can be solved explicitly. This observation
suggests that we can iteratively update one block of variables at a time while
keeping others fixed. Specifically, given the collection of
$\hat{\bm{M}}_{k}$’s, the core tensor estimate
$\hat{\mathcal{C}}=\operatorname*{arg\,min}_{\mathcal{C}}f(\mathcal{C},\\{\hat{\bm{M}}_{k}\\})$
consists of the sample averages of each tensor block. Given the block mean
$\hat{\mathcal{C}}$ and $K-1$ membership matrices, the last membership matrix
can be solved using a simple nearest neighbor search over only $R_{k}$
discrete points. The full procedure is described in Algorithm 1.

Algorithm 1 Multiway clustering based on tensor block models

1:Data tensor $\mathcal{Y}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$,
clustering size $\bm{R}=(R_{1},\ldots,R_{K})$.

2:Block mean tensor $\hat{\mathcal{C}}\in\mathbb{R}^{R_{1}\times\cdots\times
R_{K}}$, and the membership matrices $\hat{\bm{M}}_{k}$’s.

3:Initialize the marginal clustering by performing independent $k$-means on
each of the $K$ modes.

4:repeat

5: Update the core tensor
$\hat{\mathcal{C}}=\llbracket\hat{c}_{r_{1},\ldots,r_{K}}\rrbracket$.
Specifically, for each $(r_{1},\ldots,r_{K})\in[R_{1}]\times\cdots[R_{K}]$,
$\hat{c}_{r_{1},\ldots,r_{K}}={1\over
n_{r_{1},\ldots,r_{K}}}\sum_{\hat{\bm{M}}_{1}^{-1}(r_{1})\times\cdots\times\hat{\bm{M}}_{K}^{-1}(r_{K})}y_{i_{1},\ldots,i_{K}},$
(10) where $\bm{M}^{-1}_{k}(r_{k})$ denotes the indices that belong to the
$r_{k}$th cluster in the mode $k$, and
$n_{r_{1},\ldots,r_{K}}=\prod_{k}|\hat{\bm{M}}_{k}^{-1}(r_{k})|$ denotes the
number of entries in the block indexed by $(r_{1},\ldots,r_{K})$.

6: for $k$ in $\\{1,2,...,K\\}$ do

7: Update the mode-$k$ membership matrix $\hat{\bm{M}}_{k}$. Specifically, for
each $a\in[d_{k}]$, assign the cluster label $\hat{\bm{M}}_{k}(a)\in[R_{k}]$:
$\hat{\bm{M}}_{k}(a)=\operatorname*{arg\,min}_{r\in[R_{k}]}\sum_{\bm{I}_{-k}}\left(\hat{c}_{\hat{\bm{M}}_{1}(i_{1}),\ldots,r,\ldots,\hat{\bm{M}}_{K}(i_{K})}-y_{i_{1},\ldots,a,\ldots,i_{K}}\right)^{2},$
where $\bm{I}_{-k}=(i_{1},\ldots,i_{k-1},i_{k+1},\ldots,i_{K})$ denotes the
tensor coordinates except the $k$-th mode.

8: end for

9:until Convergence

Algorithm 1 can be viewed as a higher-order extension of the ordinary (one-
way) $k$-means algorithm. The core tensor $\mathcal{C}$ serves as the role of
centroids. As each iteration reduces the value of the objective function,
which is bounded below, convergence of the algorithm is guaranteed. The per-
iteration computational cost scales linearly with the sample size,
$d=\prod_{k}d_{k}$, and this complexity matches the classical tensor methods
[24, 25, 22]. We recognize that obtaining the global optimizer for such a non-
convex optimization is typically difficult [26, 1]. Following the common
practice in non-convex optimization [1], we run the algorithm multiple times,
using random initializations with independent one-way $k$-means on each of the
modes.

###  5.2 Tuning parameter selection

Algorithm 1 takes the number of clusters $\bm{R}$ as an input. In practice
such information is often unknown and $\bm{R}$ needs to be estimated from the
data $\mathcal{Y}$. We propose to select this tuning parameter using Bayesian
information criterion (BIC),

$\text{BIC}(\bm{R})=\log\left(\lVert\mathcal{Y}-\hat{\Theta}\rVert_{F}^{2}\right)+{\sum_{k}\log
d_{k}\over\prod_{k}d_{k}}p_{e},$ (11)

where $p_{e}$ is the effective number of parameters in the model. In our case
we take $p_{e}=\prod_{k}R_{k}+\sum_{k}d_{k}\log R_{k}$, which is inspired from
(7). We choose $\hat{\bm{R}}$ that minimizes $\text{BIC}(\bm{R})$ via grid
search. Our choice of BIC aims to balance between the goodness-of-fit for the
data and the degree of freedom in the population model. We test its empirical
performance in Section 7.

##  6 Extension to sparse estimation

In some large-scale applications, not every block in a data tensor is of equal
importance. For example, in the genome-wide expression data analysis, only a
few entries represent the signals while the majority come from the background
noise (see Figure 1b). While our estimator (4) is still able to handle this
scenario by assigning small values to some of the
$\hat{c}_{r_{1},\ldots,r_{K}}$’s, the estimates may suffer from high variance.
It is thus beneficial to introduce regularized estimation for better bias-
variance trade-off and improved interpretability.

Here we illustrate a sparse version of TBM by imposing regularity on block
means for localizing important blocks in the data tensor. This problem can be
formulated as a variable selection on the block parameters. We propose the
following regularized least-square estimation:

$\hat{\Theta}^{\text{sparse}}=\operatorname*{arg\,min}_{\Theta\in\mathcal{P}}\left\\{\lVert\mathcal{Y}-\Theta\rVert_{F}^{2}+\lambda\lVert\mathcal{C}\rVert_{\rho}\right\\},$

where $\mathcal{C}\in\mathbb{R}^{R_{1}\times\cdots\times R_{K}}$ is the block-
mean tensor, $\lVert\mathcal{C}\rVert_{\rho}$ is the penalty function with
$\rho$ being an index for the tensor norm, and $\lambda$ is the penalty tuning
parameter. Some widely used penalties include Lasso penalty $(\rho=1)$, sparse
subset penalty $(\rho=0)$, ridge penalty $(\rho=\text{Frobenius norm})$,
elastic net (linear combination of $\rho=1$ and $\rho=\text{Frobenius norm}$),
among many others.

For parsimony purpose, we only discuss the Lasso and sparse subset penalties;
other penalizations can be derived similarly. Sparse estimation incurs slight
changes to Algorithm 1. When updating the core tensor $\mathcal{C}$ in (10),
we fit a penalized least square problem with respect to $\mathcal{C}$. The
closed form for the entry-wise sparse estimate
$\hat{c}^{\text{sparse}}_{r_{1},\ldots,r_{K}}$ is (see Lemma 3 in the
Supplements):

$\hat{c}^{\text{sparse}}_{r_{1},\ldots,r_{K}}=\begin{cases}\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}}\mathds{1}\left\\{|\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}}|\geq{\sqrt{\lambda\over
n_{r_{1},\ldots,r_{K}}}}\right\\}&\text{if}\ \rho=0,\\\
\text{sign}(\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}})\left(|\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}}|-{\lambda\over
2n_{r_{1},\ldots,r_{K}}}\right)_{+}&\text{if}\ \rho=1,\end{cases}$ (12)

where $a_{+}=\max(a,0)$ and $\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}}$
denotes the ordinary least-square estimate in (10). The choice of penalty
$\rho$ often depends on the study goals and interpretations in specific
applications. Given a penalty function, we select the tuning parameter
$\lambda$ via BIC (11), where we modify $p_{e}$ into
$p^{\text{sparse}}_{e}=\lVert\hat{\mathcal{C}}^{\text{sparse}}\rVert_{0}+\sum_{k}d_{k}\log
R_{k}$. Here $\lVert\cdot\rVert_{0}$ denotes the number of non-zero entries in
the tensor. The empirical performance of this proposal will be evaluated in
Section 7.

##  7 Experiments

In this section, we evaluate the empirical performance of our TBM method. We
consider both non-sparse and sparse tensors, and compare the recovery accuracy
with other tensor-based methods. Unless otherwise stated, we generate Gaussian
tensors under the block model (1). The block means are generated from i.i.d.
Uniform[-3,3]. The entries in the noise tensor $\mathcal{E}$ are generated
from i.i.d. $N(0,\sigma^{2})$. In each simulation study, we report the summary
statistics across $n_{\text{sim}}=50$ replications.

###  7.1 Finite-sample performance

In the first experiment, we assess the empirical relationship between the root
mean squared error (RMSE) and the dimension. We set $\sigma=3$ and consider
tensors of order 3 and order 4 (see Figure 2). In the case of order-3 tensors,
we increase $d_{1}$ from 20 to 70, and for each choice of $d_{1}$, we set the
other two dimensions $(d_{2},d_{3})$ such that $d_{1}\log R_{1}\approx
d_{2}\log R_{2}\approx d_{3}\log R_{3}$. Recall that our theoretical analysis
suggests a convergence rate $\mathcal{O}(\sqrt{\log R_{1}/d_{2}d_{3}})$ for
our estimator. Figure 2a plots the recovery error versus the rescaled sample
size $N_{1}=\sqrt{d_{2}d_{3}/\log R_{1}}$. We find that the RMSE decreases
roughly at the rate of $1/N_{1}$. This is consistent with our theoretical
result. It is observed that tensors with a higher number of blocks tend to
yield higher recovery errors, as reflected by the upward shift of the curves
as $\bm{R}$ increases. Indeed, a higher $\bm{R}$ means a higher intrinsic
dimension of the problem, thus increasing the difficulty of the estimation.
Similar behavior can be observed in the order-4 case from Figure 2b, where the
rescaled sample size is $N_{2}=\sqrt{d_{2}d_{3}d_{4}/\log R_{1}}$.

Figure 2: Estimation error for block tensors with Gaussian noise. Each curve
corresponds to a fixed clustering size $\bm{R}$. (a) Average RMSE against
rescaled sample size $N_{1}=\sqrt{d_{2}d_{3}/\log R_{1}}$ for order-3 tensors.
(b) Average RMSE against rescaled sample size
$N_{2}=\sqrt{d_{2}d_{3}d_{4}/\log R_{1}}$ for order-4 tensors.

In the second experiment, we evaluate the selection performance of our BIC
criterion (11). Supplementary Table S1 reports the selected numbers of
clusters under various combinations of dimension $\bm{d}$, clustering size
$\bm{R}$, and noise $\sigma$. We find that, for the case $\bm{d}=(40,40,40)$
and $\bm{R}=(4,4,4)$, the BIC selection is accurate in the low-to-moderate
noise setting. In the high-noise setting with $\sigma=12$, the selected number
of clusters is slightly smaller than the true number, but the accuracy
increases when either the dimension increases to $\bm{d}=(40,40,80)$ or the
clustering size reduces to $\bm{R}=(2,3,4)$. Within a tensor, the selection
seems to be easier for shorter modes with smaller number of clusters. This
phenomenon is to be expected, since shorter mode has more effective samples
for clustering.

###  7.2 Comparison with alternative methods

Next, we compare our TBM method with two popular low-rank tensor estimation
methods: (i) CP decomposition and (ii) Tucker decomposition. Following the
literature [13, 5, 9], we perform the clustering by applying the $k$-means to
the resulting factors along each of the modes. We refer to such techniques as
CP+$k$-means and Tucker+$k$-means.

We generate noisy block tensors with five clusters on each of the modes, and
then assess both the estimation and clustering performance for each method.
Note that TBM takes a single shot to perform estimation and clustering
simultaneously, whereas CP and Tucker-based methods separate these two tasks
in two steps. We use the RMSE to assess the estimation accuracy and use the
clustering error rate (CER) to measure the clustering accuracy. The CER is
calculated using the disagreements (i.e., one minus rand index) between the
true and estimated block partitions in the three-way tensor. For fair
comparison, we provide all methods the true number of clusters.

Figure 3a shows that TBM achieves the lowest estimation error among the three
methods. The gain in accuracy is more pronounced as the noise grows. Neither
CP nor Tucker recovers the signal tensor, although Tucker appears to result in
a modest clustering performance (Figure 3b). One possible explanation is that
the Tucker model imposes orthogonality to the factors, which make the
subsequent $k$-means clustering easier than that for the CP factors. Figure
3b-c shows that the clustering error increases with noise but decreases with
dimension. This agrees with our expectation, as in tensor data analysis, a
larger dimension implies a larger sample size.

Figure 3: Performance comparison in terms of RMSE and CER. (a) Estimation
error against noise for tensors of dimension $(40,40,40)$. (b) Clustering
error against noise for tensors of dimension $(40,40,40)$. (c) Clustering
error against noise for tensors of dimension $(40,50,60)$.

Sparse case. We then evaluate the performance when the signal tensor is
sparse. The simulated model is the same as before, except that we generate
block means from a mixture of zero mass and Uniform[-3,3], with probability
$p$ (sparsity rate) and $1-p$ respectively. We generate noisy tensors of
dimension $\bm{d}=(40,40,40)$ with varying levels of sparsity and noise. We
utilize $\ell 0$-penalized TBM and primarily focus on the selection accuracy.
The performance is quantified via the the sparsity error rate, which is the
proportion of entries that were incorrectly set to zero or incorrectly set to
non-zero. We also report the proportion of true zero’s that were correctly
identified (correct zeros).

Table 2 reports the BIC-selected $\lambda$ averaged across 50 simulations. We
see a substantial benefit obtained by penalization. The proposed $\lambda$ is
able to guide the algorithm to correctly identify zero’s, while maintaining
good accuracy in identifying non-zero’s. The resulting sparsity level is close
to the ground truth. Supplementary Figure S1 shows the estimation error and
sparsity error against $\sigma$ when $p=0.8$. Again, the sparse TBM
outperforms the other methods.

Sparsity ($p$) | Noise $(\sigma)$ | BIC-selected $\lambda$ | Estimated Sparsity Rate | Correct Zero Rate | Sparsity Error Rate
---|---|---|---|---|---
0.5 | 4 | 136.0(37.5) | 0.55(0.04) | 1.00(0.02) | 0.06(0.03)
0.5 | 8 | 439.2(80.2) | 0.58(0.06) | 0.94(0.08) | 0.15(0.07)
0.8 | 8 | 458.0(63.3) | 0.81(0.15) | 0.87(0.16) | 0.21(0.13)

Table 2: Sparse TBM for estimating tensors of dimension $\bm{d}=(40,40,40)$.
The reported statistics are averaged across 50 simulations with standard
deviation given in parentheses. Number in bold indicates the ground truth is
within 2 standard deviations of the sample average.

###  7.3 Real data analysis

Lastly, we apply our method on two real datasets. We briefly summarize the
main findings here; the detailed information can be found in the Supplements.

The first dataset is a real-valued tensor, consisting of approximate 1 million
expression values from 13 brain tissues, 193 individuals, and 362 genes [4].
We subtracted the overall mean expression from the data, and applied the $\ell
0$-penalized TBM to identify important blocks in the resulting tensor. The top
blocks exhibit a clear tissues $\times$ genes specificity (Supplementary Table
S2). In particular, the top over-expressed block is driven by tissues
{_Substantia nigra, Spinal cord_} and genes {_GFAP, MBP_}, suggesting their
elevated expression across individuals. In fact, _GFAP_ encodes filament
proteins for mature astrocytes and _MBP_ encodes myelin sheath for
oligodendrocytes, both of which play important roles in the central nervous
system [27]. Our method also identifies blocks with extremely negative means
(i.e. under-expressed blocks). The top under-expressed block is driven by
tissues {_Cerebellum, Cerebellar Hemisphere_} and genes {_CDH9, GPR6, RXFP1,
CRH, DLX5/6, NKX2-1, SLC17A8_}. The gene _DLX6_ encodes proteins in the
forebrain development [27], whereas cerebellum tissues are located in the
hindbrain brain. The opposite spatial function is consistent with the observed
under-expression pattern.

The second dataset we consider is the _Nations_ data [2]. This is a $14\times
14\times 56$ binary tensor consisting of 56 political relationships of 14
countries between 1950 and 1965. We note that 78.9% of the entries are zero.
Again, we applied the $\ell 0$-penalized TBM to identify important blocks in
the data. We found that the 14 countries are naturally partitioned into 5
clusters, two representing neutral countries {_Brazil, Egypt, India, Israel,
Netherlands_} and {_Burma, Indonesia, Jordan_}, one eastern bloc {_China,
Cuba, Poland, USSA_}, and two western blocs, {_USA_} and {_UK_}. The relation
types are partitioned into 7 clusters, among which the exports-related
activities {_reltreaties, book translations, relbooktranslations, exports3,
relexporsts_} and NGO-related activities {_relintergovorgs, relngo,
intergovorgs3, ngoorgs3_} are two major clusters that involve the connection
between neutral and western blocs. Other top blocks are described in the
Supplement.

We compared the goodness-of-fit of various clustering methods on the _Brain
expression_ and _Nations_ datasets. Because the code of CoCo method [13] is
not yet available, we excluded it from our numerical comparison (See Section 4
for the theoretical comparison with CoCo). Table 3 summarizes the proportion
of variance explained by each clustering method:

Dataset | TBM | TBM-sparse | CP+$k$-means | Tucker+$k$-means | CoTeC [12]
---|---|---|---|---|---
Brain expression | 0.856 | 0.855 | 0.576 | 0.434 | 0.849
Nations | 0.439 | 0.433 | 0.324 | 0.253 | 0.419

Table 3: Comparison of goodness-of-fit in the _Brain_ expression and _Nations_
datasets.

Our method (TBM) achieves the highest variance proportion, suggesting that the
entries within the same cluster are close (i.e., a good clustering). As
expected, the sparse TBM results in a slightly lower proportion, because it
has a lower model complexity at the cost of small bias. It is remarkable that
the sparse TBM still achieves a higher goodness-of-fit than others. The
improved interpretability with little loss of accuracy makes the sparse TBM
appealing in applications.

##  8 Conclusion

We have developed a statistical setting for studying the tensor block model.
Under the assumption that tensor entries are distributed with a block-specific
mean, our estimator achieves a convergence rate $\mathcal{O}(\sum_{k}d_{k}\log
R_{k})$ which is faster than previously possible. Our TBM method applies to a
broad range of data distributions and can handle both sparse and dense data
tensor. We demonstrate the benefit of sparse regularity in power of detection.
In specific applications, prior knowledge may suggest other regularities for
parameters. For example, in the multi-layer network analysis, sometimes it may
be reasonable to impose symmetry on the parameters along certain modes. In
some other applications, non-negativity of parameter values may be enforced.
We leave these directions for future study.

## Acknowledgements

This research was supported by NSF grant DMS-1915978 and the University of
Wisconsin-Madison, Office of the Vice Chancellor for Research and Graduate
Education with funding from the Wisconsin Alumni Research Foundation.

##  Appendix A Proofs

###  A.1 Stochastic tensor block model

The following property shows that Bernoulli distribution belongs to the sub-
Gaussian family with a subgaussianity parameter $\sigma$ equal to $1/4$.

######  Property 1.

Suppose $x\sim\text{Bernoulli}(p)$, then $x\sim\text{sub-Gaussian}({1\over
4})$.

###### Proof.

For all $\lambda\in\mathbb{R}$, we have

$\ln(\mathbb{E}(e^{\lambda(x-p)})=\ln\left(pe^{\lambda(1-p)}+(1-p)e^{-p\lambda}\right)=-p\lambda+\ln(1+pe^{\lambda}-p)\leq{\lambda^{2}\over
8}.$

Therefore $\mathbb{E}(e^{\lambda(x-p)})\leq e^{\lambda^{2}(1/4)/2}$. ∎

###  A.2 Proof of Proposition 1

###### Proof.

Let $\mathbb{P}_{\Theta}$ denotes the (either Gaussian or Bernoulli) tensor
block model, where
$\Theta=\mathcal{C}\times_{1}\bm{M}_{1}\times_{2}\cdots\times_{K}\bm{M}_{K}$
parameterizes the mean tensor. Since the mapping
$\Theta\mapsto\mathbb{P}_{\Theta}$ is one-to-one, $\Theta$ is identifiable.
Now suppose that $\Theta$ can be decomposed in two ways,
$\Theta=\Theta(\\{\bm{M}_{k}\\},\mathcal{C})=\Theta(\\{\tilde{\bm{M}}_{k}\\},\tilde{\mathcal{C}})$.
Based on the Assumption 1, we have

$\Theta=\mathcal{C}\times_{1}\bm{M}_{1}\times_{2}\cdots\times_{K}\bm{M}_{K}=\tilde{\mathcal{C}}\times_{1}\tilde{\bm{M}}_{1}\times_{2}\cdots\times_{K}\tilde{\bm{M}}_{K},$
(13)

where $\mathcal{C}$,
$\tilde{\mathcal{C}}\in\mathbb{R}^{R_{1}\times\cdots\times R_{K}}$ are two
irreducible cores, and $\bm{M}_{k},\tilde{\bm{M}}_{k}\in\\{0,1\\}^{R_{k}\times
d_{k}}$ are membership matrices for all $k\in[K]$. We will prove by
contradiction that $\bm{M}_{k}$ and $\tilde{\bm{M}}_{k}$ induce the same
partition of $[d_{k}]$, for all $k\in[K]$.

Suppose the above claim does not hold. Then there exists a mode $k\in[K]$ such
that the $\bm{M}_{k},\tilde{\bm{M}}_{k}$ induce two different partitions of
$[d_{k}]$. Without loss of generality, we assume $k=1$. The definition of
partition implies that there exists a pair of indices $i\neq j$,
$i,j\in[d_{1}]$, such that, $i,j$ belong to the same cluster based on
$\bm{M}_{1}$, but they belong to different clusters based on
$\tilde{\bm{M}}_{1}$. Let
$\mathcal{A}\neq\mathcal{B},\mathcal{A},\mathcal{B}\subset[d_{1}]$
respectively denote the clusters that $i$ and $j$ belong to, based on
$\tilde{\bm{M}}_{1}$. The left-hand side of (13) implies

$\Theta_{i,i_{2},\ldots,i_{K}}=\Theta_{j,i_{2},\ldots,i_{K}},\quad\text{for
all }(i_{2},\ldots,i_{K})\in[d_{2}]\times\cdots\times[d_{K}].$ (14)

On the other hand, (13) implies

$\Theta_{i,i_{2},\ldots,i_{K}}=\Theta_{k,i_{2},\ldots,i_{K}},\quad\text{for
all }k\in\mathcal{A}\text{ and all
}(i_{2},\ldots,i_{K})\in[d_{2}]\times\cdots\times[d_{K}],$ (15)

and

$\Theta_{j,i_{2},\ldots,i_{K}}=\Theta_{k,i_{2},\ldots,i_{K}},\quad\text{for
all }k\in\mathcal{B}\text{ and all
}(i_{2},\ldots,i_{K})\in[d_{2}]\times\cdots\times[d_{K}].$ (16)

Combining (14), (15) and (16), we have

$\Theta_{i,i_{2},\ldots,i_{K}}=\Theta_{k,i_{2},\ldots,i_{K}},\quad\text{for
all }k\in\mathcal{A}\cup\mathcal{B}\text{ and all
}(i_{2},\ldots,i_{K})\in[d_{2}]\times\cdots\times[d_{K}].$ (17)

Equation (17) implies that $\mathcal{A}$ and $\mathcal{B}$ can be merged into
one cluster. This contradicts the irreducibility assumption of the core tensor
$\tilde{\mathcal{C}}$. Therefore, $\bm{M}_{1}$ and $\tilde{\bm{M}}_{1}$ induce
a same partition of $[d_{1}]$, and thus they are equal up to permutation of
cluster labels. The proof is now complete. ∎

###  A.3 Proof of Theorem 1

The following lemma is useful for the proof of Theorem 1.

######  Lemma 1.

Suppose $\mathcal{Y}=\Theta_{\text{true}}+\mathcal{E}$ with
$\Theta_{\text{true}}\in\mathcal{P}$. Let
$\hat{\Theta}=\arg\min_{\Theta\in\mathcal{P}}\lVert\hat{\Theta}-\mathcal{Y}\rVert_{F}^{2}$
be the least-square estimator of $\Theta_{\text{true}}$. We have

$\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}\leq
2\sup_{\mu\in{\mathcal{P}-\mathcal{P}^{\prime}\over|\mathcal{P}-\mathcal{P}^{\prime}|}}\langle\mu,\mathcal{E}\rangle,$

where
$\mathcal{P}-\mathcal{P}^{\prime}=\\{\Theta-\Theta^{\prime}\colon\Theta,\Theta^{\prime}\in\mathcal{P}\\}$
and $\mathcal{S}/|\mathcal{S}|=\\{s/\lVert s\rVert_{2}\colon
s\in\mathcal{S}\\}$.

###### Proof.

Based on the definition of least-square estimator, we have

$\lVert\hat{\Theta}-\mathcal{Y}\rVert_{F}^{2}\leq\lVert\Theta_{\text{true}}-\mathcal{Y}\rVert_{F}^{2}.$
(18)

Combining (18) with the fact

$\begin{split}\lVert\hat{\Theta}-\mathcal{Y}\rVert_{F}^{2}&=\lVert\hat{\Theta}-\Theta_{\text{true}}+\Theta_{\text{true}}-\mathcal{Y}\rVert_{F}^{2}\\\
&=\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}^{2}+\lVert\Theta_{\text{true}}-\mathcal{Y}\rVert_{F}^{2}+2\langle\hat{\Theta}-\Theta_{\text{true}},\Theta_{\text{true}}-\mathcal{Y}\rangle,\end{split}$
(19)

yields

$\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}^{2}\leq
2\langle\hat{\Theta}-\Theta_{\text{true}},\mathcal{Y}-\Theta_{\text{true}}\rangle=2\langle\hat{\Theta}-\Theta_{\text{true}},\mathcal{E}\rangle.$

Dividing each side by $\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}$, we
have

$\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}\leq
2\left\langle{\hat{\Theta}-\Theta_{\text{true}}\over\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}},\mathcal{E}\right\rangle.$

The desired inequality follows by noting
${\hat{\Theta}-\Theta_{\text{true}}\over\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}}\in{\mathcal{P}-\mathcal{P}^{\prime}\over|\mathcal{P}-\mathcal{P}^{\prime}|}$.
∎

###### Proof of Theorem 1.

To study the performance of the least-square estimator $\hat{\Theta}$, we need
to introduce some additional notation. We view the membership matrix
$\bm{M}_{k}$ as an onto function $\bm{M}_{k}\colon[d_{k}]\mapsto[R_{k}]$. With
a little abuse of notation, we still use $\bm{M}_{k}$ to denote the mapping
function and write $\bm{M}_{k}\in R_{k}^{d_{k}}$ by convention. We use
$\bm{M}=\\{\bm{M}_{k}\\}_{k\in[K]}$ to denote the collection of $K$ membership
matrices, and write $\mathcal{M}=\\{\bm{M}\colon\text{$\bm{M}$ is the
collection of membership matrices $\bm{M}_{k}$'s}\\}$. For any set $J$, $|J|$
denotes its cardinality. Note that $|\mathcal{M}|\leq\prod_{k}R_{k}^{d_{k}}$,
because each $\bm{M}_{k}$ can be identified by a partition of $[d_{k}]$ into
$R_{k}$ disjoint non-empty sets.

For ease of notation, we define $d=\prod_{k}d_{k}$ and $R=\prod_{k}R_{k}$. We
sometimes identify a tensor in $\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$
with a vector in $\mathbb{R}^{d}$. By the definition of the parameter space
$\mathcal{P}$, the element $\Theta\in\mathcal{P}$ can be equivalently
identified by $\Theta=\Theta(\bm{M},\bm{C})$, where $\bm{M}\in\mathcal{M}$ is
the collection of $K$ membership matrices and
$\bm{C}=\text{vec}(\mathcal{C})\in\mathbb{R}^{R}$ is the core tensor. Note
that, for a fixed clustering structure $\bm{M}$, the space consisting of
$\Theta=\Theta(\bm{M},\cdot)$ is a linear space of dimension $R$.

Now consider the least-square estimator

$\hat{\Theta}=\operatorname*{arg\,min}_{\Theta\in\mathcal{P}}\\{-2\langle\mathcal{Y},\Theta\rangle+\lVert\Theta\rVert_{F}^{2}\\}=\operatorname*{arg\,min}_{\Theta\in\mathcal{P}}\\{\lVert\mathcal{Y}-\Theta\rVert_{F}^{2}\\}.$

Based on the Lemma 1,

$\begin{split}\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}&\leq
2\sup_{\Theta\in\mathcal{P}}\sup_{\Theta^{\prime}\in\mathcal{P}}\Big{\langle}{\Theta-\Theta^{\prime}\over\lVert\Theta-\Theta^{\prime}\rVert_{F}},\mathcal{E}\Big{\rangle}\\\
&\leq
2\sup_{\bm{M},\bm{M}^{\prime}\in\mathcal{M}}\sup_{\bm{C},\bm{C}^{\prime}\in\mathbb{R}^{R}}\Big{\langle}{\Theta(\bm{M},\bm{C})-\Theta^{\prime}(\bm{M}^{\prime},\bm{C}^{\prime})\over\lVert\Theta(\bm{M},\bm{C})-\Theta^{\prime}(\bm{M}^{\prime},\bm{C}^{\prime})\rVert_{F}},\mathcal{E}\Big{\rangle}.\end{split}$
(20)

By union bound, we have, for any $t>0$,

$\begin{split}\mathbb{P}\left(\lVert\hat{\Theta}-\Theta_{\text{true}}\rVert_{F}>t\right)&\leq\mathbb{P}\left(\sup_{\bm{M},\bm{M}^{\prime}\in\mathcal{M}}\sup_{\bm{C},\bm{C}^{\prime}\in\mathbb{R}^{R}}\left|\Big{\langle}{\Theta(\bm{M},\bm{C})-\Theta^{\prime}(\bm{M}^{\prime},\bm{C}^{\prime})\over\lVert\Theta(\bm{M},\bm{C})-\Theta^{\prime}(\bm{M}^{\prime},\bm{C}^{\prime})\rVert_{F}},\
\mathcal{E}\Big{\rangle}\right|>{t\over 2}\right)\\\
&\leq\sum_{\bm{M},\bm{M}^{\prime}\in\mathcal{M}}\mathbb{P}\left(\sup_{\bm{C}^{\prime}\in\mathbb{R}^{R}}\sup_{\bm{C}\in\mathbb{R}^{R}}\left|\Big{\langle}{\Theta(\bm{M},\mathcal{C})-\Theta^{\prime}(\bm{M}^{\prime},\mathcal{C})\over\lVert\Theta(\bm{M},\mathcal{C})-\Theta^{\prime}(\bm{M}^{\prime},\mathcal{C})\rVert_{F}},\
\mathcal{E}\Big{\rangle}\right|\geq{t\over 2}\right)\\\
&\leq|\mathcal{M}|^{2}C_{1}^{R}\exp\left(-{C_{2}t^{2}\over
32\sigma^{2}}\right)\\\ &=\exp\left(2\sum_{k}d_{k}\log
R_{k}+C_{1}\prod_{k}R_{k}-{C_{2}t^{2}\over 32\sigma^{2}}\right),\end{split}$
(21)

for two universal constants $C_{1},C_{2}>0$. Here the third line follows from
[28] (Theorem 1.19) and the fact that $\Theta=\Theta(\bm{M},\cdot)$ lies in a
linear space of dimension $R$. The last line uses
$|\mathcal{M}|\leq\prod_{k}R_{k}^{d_{k}}$ and $R=\prod_{k}R_{k}$. Choosing
$t=C\sigma\sqrt{\prod_{k}R_{k}+\sum_{k}d_{k}\log R_{k}}$ yields the desired
bound. ∎

###  A.4 Proof of Theorem 2

First we give a list of notation used in the proof. For ease of notation, we
allow the basic arithmetic operators ($+,-,\geq$, etc) to be applied to pairs
of vectors in an element-wise manner.

####  A.4.1 Notations

$\bm{M}_{k}=\llbracket m_{ir}^{(k)}\rrbracket\in\\{0,1\\}^{d_{k}\times
R_{k}}$: the mode-$k$ membership matrix. The element $m_{ir}^{(k)}=1$ if and
only if the $i$th slide in mode $k$ belongs to the $r$th cluster.

$\bm{M}_{k,\text{true}}$, $\hat{\bm{M}}_{k}\in\\{0,1\\}^{d_{k}\times R_{k}}$:
the true and estimated mode-$k$ cluster membership matrices, respectively.

$\bm{p}^{(k)}=\llbracket p^{(k)}_{r}\rrbracket\in[0,1]^{R_{k}}$: the marginal
cluster proportion vector listing the relative cluster sizes along the mode
$k$. The element
$p_{r}^{(k)}=\frac{1}{d_{k}}\sum_{i=1}^{d_{k}}\mathds{1}\\{m_{ir}^{(k)}=1\\}$
denotes the proportion of the $r$th cluster. The cluster proportion vector
$\bm{p}^{(k)}=\bm{p}^{(k)}(\bm{M}_{k})$ can be viewed as a function of
$\bm{M}_{k}$.

$\bm{p}^{(k)}_{\text{true}}$, $\hat{\bm{p}}^{(k)}\in[0,1]^{R_{k}}$: the true
and estimated mode-$k$ cluster proportion vectors, respectively.

$\bm{D}^{(k)}=\llbracket D^{(k)}_{rr^{\prime}}\rrbracket\in[0,1]^{R_{k}\times
R_{k}}$: the mode-$k$ confusion matrix between clustering
$\bm{M}_{k,\text{true}}$ and $\hat{\bm{M}}_{k}$. The entries in the confusion
matrix is
$D_{rr^{\prime}}^{(k)}=\frac{1}{d_{k}}\sum_{i=1}^{d_{k}}\mathbb{I}\\{m_{ir,\text{true}}^{(k)}=\hat{m}_{ir^{\prime}}^{(k)}=1\\}$.
The confusion matrix $\bm{D}^{(k)}={1\over
d_{k}}\bm{M}^{T}_{k,\text{true}}\hat{\bm{M}}_{k}$ is a function of
$\bm{M}_{k,\text{true}}$ and $\hat{\bm{M}}_{k}$.

$\mathcal{J}_{\tau}=\\{(\bm{M}_{1},\ldots,\bm{M}_{K}):\bm{p}^{(k)}(\bm{M}_{k})\geq\tau\text{
for all }k\in[K]\\}$: the set of all possible partitions that satisfy the
marginal non-degenerating assumption.

$\mathcal{I}\subset 2^{[d_{1}]}\times\cdot\cdot\cdot\times 2^{[d_{K}]}$: the
set of blocks that satisfy the marginal non-degenerating assumption for all
$k\in[K]$;

$L=\inf\\{|I|:I\in\mathcal{I}\\}$: the minimum block size in $\mathcal{I}$.

$\lVert\mathcal{A}\rVert_{\max}=\max_{r_{1},\ldots,r_{K}}|a_{r_{1},\ldots,r_{K}}|$
for any tensor $\mathcal{A}=\llbracket
a_{i_{1},\ldots,i_{K}}\rrbracket\in\mathbb{R}^{R_{1}\times\ldots\times
R_{K}}$.

$f(x)=x^{2}$: the quadratic objective function.

######  Remark 1.

By definition, the confusion matrix $\bm{D}^{(k)}$ satisfies the following two
properties:

  1. 1.

$\bm{D}^{(k)}\bm{1}=\bm{p}^{(k)}_{\text{true}}$,
$(\bm{D}^{(k)})^{T}\bm{1}=\hat{\bm{p}}^{(k)}$.

  2. 2.

The estimated clustering matches the true clustering if and only if
$\bm{D}^{(k)}$ equals to the diagonal matrix up to permutation.

####  A.4.2 Auxiliary Results

Recall that the objective function in our tensor block model is

$\displaystyle f(\mathcal{C},\\{\bm{M}_{k}\\})=\langle\mathcal{Y},\
\Theta\rangle-\frac{\lVert\Theta\rVert_{F}^{2}}{2},$ (22)
$\displaystyle\text{where
}\Theta=\mathcal{C}\times_{1}\bm{M}_{1}\times_{2}\cdots\times_{K}\bm{M}_{K},$
(23)

where $\mathcal{Y}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$ is the data,
$\mathcal{C}$ is the core tensor of interest, and $\\{\bm{M}_{k}\\}$ is the
membership matrices of interest. Without loss of generality, we will work with
the scaled objective ${2\over\prod_{k}d_{k}}f(\mathcal{C},\\{\bm{M}_{k}\\})$.
With a little abuse of notation, we still denote the scaled function as
$f(C,\\{\bm{M}_{k}\\})$.

We will prove that, if there is non-negligible mismatch between
$\\{\hat{\bm{M}}_{k}\\}$ and $\\{\bm{M}_{k,\text{true}}\\}$, then
$\\{\hat{\bm{M}}_{k}\\}$ cannot be the optimizer to (22). To show this, we
investigate the objective values at the global optimizer vs. at the true
parameter. The deviation between these two values comes from two aspects: the
label assignments (i.e., the estimation of $\\{\bm{M}_{k}\\}$) and the
estimation of the core tensor. In what follows, we tease apart these two
aspects.

  1. 1.

First, suppose the partitions $\\{\bm{M}_{k}\\}$ are given, which are not
necessarily equal to $\\{\bm{M}_{k,\text{true}}\\}$. We now assess the
stochastic error due to estimation of $\mathcal{C}$, conditional on
$\\{\bm{M}_{k}\\}$. In such a case, the core
$\hat{\mathcal{C}}=\operatorname*{arg\,min}_{\mathcal{C}}f(\mathcal{C},\\{\bm{M}_{k}\\})$
can be solved explicitly. Specifically, the optimizer
$\hat{\mathcal{C}}=\llbracket\hat{c}_{r_{1},\ldots,r_{K}}\rrbracket$ consists
of the sample averages of each tensor block, where

$\displaystyle\displaystyle\hat{c}_{r_{1},\ldots,r_{K}}$
$\displaystyle=\hat{c}_{r_{1},\ldots,r_{K}}(\\{\bm{M}_{k}\\})$ (24)
$\displaystyle={1\over d_{1}\cdots d_{K}}\frac{1}{p_{r_{1}}^{(1)}\cdots
p_{r_{K}}^{(K)}}\left[\mathcal{Y}\times_{1}\bm{M}^{T}_{1}\times_{2}\cdots\times_{K}\bm{M}^{T}_{K}\right]_{r_{1},\ldots,r_{K}}$
(25)

where the marginal cluster proportion $p^{(k)}_{r_{k}}$ is induced by the
clustering $\bm{M}_{k}$.

Define a new cost function
$F(\bm{M}_{1},\ldots,\bm{M}_{K})=-f(\hat{\mathcal{C}},\bm{M}_{1},\ldots,\bm{M}_{k})$,
where $\hat{\mathcal{C}}=\llbracket\hat{c}_{i_{1},\ldots,i_{K}}\rrbracket$ is
expressed in (24). A straightforward calculation shows that the function
$F(\cdot)$ has the form

$F(\bm{M}_{1},\ldots,\bm{M}_{K})=\sum_{r_{1},\ldots,r_{K}}\left(\prod_{k}p^{(k)}_{r_{k}}\right)\hat{c}_{r_{1},\ldots,r_{K}}^{2}.$
(26)

Let
$G(\bm{M}_{1},\ldots,\bm{M}_{k})=\mathbb{E}(F(\bm{M}_{1},\ldots,\bm{M}_{K}))$,
where the expectation is taken with respect to the
$\hat{\mathcal{C}}=\llbracket\hat{c}_{r_{1},\ldots,r_{K}}\rrbracket$. We have
that

$G(\bm{M}_{1},\ldots,\bm{M}_{K})=\sum_{r_{1},\ldots,r_{K}}\left(\prod_{k}p^{(k)}_{r_{k}}\right)\mu^{2}_{r_{1},\ldots,r_{K}},$
(27)

where

$\mu_{r_{1},\ldots,r_{K}}=\mathbb{E}(\hat{c}_{r_{1},\ldots,r_{K}})={1\over\prod_{k}p^{(k)}_{r_{k}}}\left[\mathcal{C}\times_{1}\bm{D}^{(1)^{T}}\times_{2}\cdots\times_{K}\bm{D}^{(K)^{T}}\right]_{r_{1},\ldots,r_{K}}$

is the expectation of the average of $y_{i_{1},\ldots,i_{K}}$ over the tensor
block indexed by $(r_{1},\ldots,r_{K})$, and $\bm{D}^{(k)}=\llbracket
D^{(k)}_{i_{k}j_{k}}\rrbracket$ is the confusion matrix between
$\bm{M}_{k,\text{true}}$ and $\bm{M}_{k}$.

The deviation
$F(\bm{M}_{1},\ldots,\bm{M}_{K})-G(\bm{M}_{1},\ldots,\bm{M}_{K})$ quantifies
the stochastic error caused by the core tensor estimation. We sometimes use
$G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})$ to denote
$G(\bm{M}_{1},\ldots,\bm{M}_{k})$ if we want to emphasize the error caused by
mismatch in label assignments. Based on (26) and (27), we define a residual
tensor for the block means:

$\displaystyle\mathcal{R}(\bm{M}_{1},\ldots,\bm{M}_{K})=\llbracket
R_{r_{1},\ldots,r_{K}}\rrbracket,\text{where}$ (28) $\displaystyle
R_{r_{1},\ldots,r_{K}}=\hat{c}_{r_{1},\ldots,r_{K}}-\mu_{r_{1},\ldots,r_{K}},\quad\text{for
all }(r_{1},\ldots,r_{K})\in[R_{1}]\times\cdots\times[R_{K}].$ (29)

Note that, conditional on $\\{\bm{M}_{k}\\}$, the entries
$R_{r_{1},\ldots,r_{K}}$ in the residual tensor are independent sub-Gaussian
with parameter depending on the size of the $(r_{1},\ldots,r_{K})$th block.

  2. 2.

Next, we free $\\{\bm{M}_{k}\\}$ and quantify the total stochastic deviation.
Note that optimizing (22) is equivalent to optimizing (26) with respect to
$\\{\bm{M}_{k}\\}$. So the least-square estimator of $\\{\bm{M}_{k}\\}$ can be
expressed as

$(\hat{\bm{M}}_{1},\ldots,\hat{\bm{M}}_{K})=\displaystyle\operatorname*{arg\,max}_{(\bm{M}_{1},\ldots,\bm{M}_{K})\in\mathcal{J}_{\tau}}F(\bm{M}_{1},\ldots,\bm{M}_{K}).$
(30)

The expectation (with respect to $\hat{\mathcal{C}}$) of the objective value
at the true parameter is

$G(\bm{M}_{1,\text{true}},\
\ldots,\bm{M}_{K,\text{true}})=\sum_{r_{1},\ldots,r_{K}}p^{(1)}_{r_{1},\text{true}}\cdots
p^{(K)}_{r_{K},\text{true}}c^{2}_{r_{1},\ldots,r_{K},\text{true}}$

We use $G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})-G(\bm{M}_{1,\text{true}},\
\ldots,\bm{M}_{K,\text{true}})$ to measure the stochastic deviation caused by
mismatch in label assignments; and use
$F(\bm{M}_{1},\ldots,\bm{M}_{K})-G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})$ to
measure stochastic deviation caused by estimation of core tensors.

The following lemma shows that, if there is non-negligible mismatch between
$\bm{M}_{k,\text{true}}$ and $\hat{\bm{M}}_{k}$, then $\hat{\bm{M}}_{k}$
cannot be the global optimizer to the objective function (22).

######  Lemma 2.

Consider partitions that satisfying
$(\bm{M}_{1},\ldots,\bm{M}_{K})\in\mathcal{J}_{\tau}$, for some $\tau>0$.
Define the minimal gap between block means $\delta^{(k)}=\min_{r_{k}\neq
r_{k}^{\prime}}$
$\max_{r_{1},\ldots,r_{k-1},r_{k+1},\ldots,r_{K}}(c_{r_{1},\ldots,r_{k},\ldots,r_{K}}-c_{r_{1},\ldots,r_{k}^{\prime},\ldots,r_{K}})^{2}>0$
and assume $\delta_{\min}=\min_{k}\delta^{(k)}>0$. For any fixed
$\varepsilon>0$, suppose
$\text{MCR}(\bm{M}_{k,\text{true}},\hat{\bm{M}}_{k})\geq\varepsilon$ for some
$k\in[K]$. Then, we have

$G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})-G(\bm{M}_{1,\text{true}},\
\ldots,\bm{M}_{K,\text{true}})\leq-\frac{1}{4}\varepsilon\tau^{K-1}\delta_{\min},$

where $\bm{D}^{(k)}$ is the confusion matrix between $\bm{M}_{k,\text{true}}$
and $\hat{\bm{M}}_{k}$.

###### Proof of Lemma 2.

For ease of notation, we drop the subscript “true” and simply write
$p^{(k)}_{r_{k}}$, $\bm{M}_{k}$, $\mathcal{C}$, etc. as the true parameters.
The corresponding estimators are denoted as $\hat{p}^{(k)}_{r_{k}}$,
$\hat{\bm{M}}_{k}$, etc. Recall that

$G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})=\sum_{r_{1},\ldots,r_{K}}\hat{p}^{(1)}_{r_{1}}\cdots\hat{p}^{(K)}_{r_{K}}\mu^{2}_{r_{1},\ldots,r_{K}},$
(31)

where $\hat{p}^{(k)}_{r_{k}}$ is the marginal cluster proportion induced by
$\hat{\bm{M}}_{k}$, and $\mu_{r_{1},\ldots,r_{K}}$ is the expected block mean
induced by $\hat{\bm{M}}_{k}$:

$\mu_{r_{1},\ldots,r_{K}}=\mu_{r_{1},\ldots,r_{K}}(\hat{\bm{M}}_{1},\ldots,\hat{\bm{M}}_{K})={1\over\prod_{k}\hat{p}^{(k)}_{r_{k}}}\left[\mathcal{C}\times_{1}\bm{D}^{(1)^{T}}\times_{2}\cdots\times_{K}\bm{D}^{(K)^{T}}\right]_{r_{1},\ldots,r_{K}}.$

We provide the proof for $k=1$. The proof for other $k\in[K]$ is similar. The
condition on MCR implies that, there exist some $r_{1}\in[R_{1}]$ and some
$a_{1}\neq a_{1}^{\prime}\in[R_{1}]$, such that
$\min\\{D_{a_{1}r_{1}}^{(1)},D_{a_{1}^{\prime}r_{1}}^{(1)}\\}\geq\varepsilon$.
Because the minimal gap between tensor block means are non-zero, we choose
$(a_{2},\ldots,a_{K})$ such that
$(c_{a_{1},a_{2},\ldots,a_{K}}-c_{a_{1}^{\prime},a_{2},\ldots,a_{K}})^{2}=\displaystyle\max_{a_{2},\ldots,a_{K}}(c_{a_{1},a_{2},\ldots,a_{K}}-c_{a_{1}^{\prime},a_{2},\ldots,a_{K}})^{2}>0$.

Let $\mathcal{N}=\llbracket
c^{2}_{a_{1},\ldots,a_{K}}\rrbracket\in\mathbb{R}^{R_{1}\times\cdots\times
R_{K}}$ be the quadratic loss evaluated at block,
$W_{r_{1},\ldots,r_{K}}=\prod_{k}\hat{p}^{(k)}_{r_{k}}>0$ the size for the
block indexed by $(r_{1},\ldots,r_{K})$. For ease of notation, we drop the
subscript ${(r_{1},\ldots,r_{K})}$ and simply write $W$.

Based on the convexity of quadratic loss, there exists $c_{*}\in\mathbb{R}$
such that the weighted quadratic loss can be expressed as

$\displaystyle[\mathcal{N}\times_{1}\bm{D}^{(1)^{T}}\times_{2}\cdots\times_{K}\bm{D}^{(K)^{T}}]_{r_{1},\ldots,r_{K}}$
(32) $\displaystyle=$ $\displaystyle\
D_{a_{1}r_{1}}^{(1)}D_{a_{2}r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}c^{2}_{a_{1},a_{2},\ldots.,a_{K}}+D_{a_{1}^{\prime}r_{1}}^{(1)}D_{a_{2}r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}c^{2}_{a_{1}^{\prime},a_{2},\ldots,a_{K}}+$ (33)
$\displaystyle(W-D_{a_{1}r_{1}}^{(1)}D_{a_{2}r_{2}}^{(1)}\cdots
D_{a_{K}r_{K}}^{(K)}-D_{a_{1}^{\prime}r_{1}}^{(1)}D_{a_{2}r_{2}}^{(1)}\cdots
D_{a_{K}r_{K}}^{(K)})c^{2}_{*}.$ (34)

Recall that $\mu_{r_{1},\ldots,r_{K}}={1\over
W}[\mathcal{C}\times_{1}\bm{D}^{(1)^{T}}\times_{2}\cdots\times_{K}\bm{D}^{(K)^{T}}]_{r_{1},\ldots,r_{K}}$
is the $(r_{1},\ldots,r_{k})$-th weighted entry of the block means. By the
Taylor expansion of quadratic loss function at $\mu_{r_{1},\ldots,r_{K}}$, we
have

$\displaystyle{1\over
W}[\mathcal{N}\times_{1}\bm{D}^{(1)^{T}}\times_{2}\cdots\times_{K}\bm{D}^{(K)^{T}}]_{r_{1},\ldots,r_{K}}-\mu^{2}_{r_{1},\ldots,r_{K}}$
$\displaystyle\geq$ $\displaystyle\ {1\over
2W}D_{a_{1}r_{1}}^{(1)}D_{a_{2}r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}(c_{a_{1},a_{2},\ldots,a_{K}}-\mu_{r_{1},\ldots,r_{K}})^{2}+$
$\displaystyle{1\over
2W}D_{a^{\prime}_{1},r_{1}}^{(1)}D_{a_{2}r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}(c_{a^{\prime}_{1},a_{2},\ldots,a_{K}}-\mu_{r_{1},\ldots,r_{K}})^{2}+$
$\displaystyle\frac{1}{2W}\left(W-D_{a_{1}r_{1}}^{(1)}D_{a_{2},r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}-D_{a^{\prime}_{1},r_{1}}^{(1)}D_{a_{2},r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}\right)(c_{*}-\mu_{r_{1},\ldots,r_{K}})^{2}.$ (35)

Combining (A.4.2) and basic inequality $(a^{2}+b^{2})\geq{1\over 2}(a+b)^{2}$
gives

$\displaystyle\ {1\over
W}[\mathcal{N}\times_{1}\bm{D}^{(1)^{T}}\times_{2}\cdots\times_{K}\bm{D}^{(K)^{T}}]_{r_{1},\ldots,r_{K}}-\mu_{r_{1},\ldots,r_{K}}^{2}$
$\displaystyle\geq$
$\displaystyle\frac{1}{4W}\min\left\\{D_{a_{1}r_{1}}^{(1)},\
D_{a_{1}^{\prime}r_{1}}^{(1)}\right\\}D_{a_{2}r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}(c_{a_{1},\ldots,a_{K}}-c_{a_{1}^{\prime},\ldots,a_{K}})^{2}$
$\displaystyle\geq$ $\displaystyle\ \frac{\varepsilon
D_{a_{2}r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}}{4W}(c_{a_{1},a_{2},\ldots,a_{K}}-c_{a_{1}^{\prime},a_{2},\ldots,a_{K}})^{2}.$
(36)

The inequality (A.4.2) only holds for a certain $r_{1}\in[R_{1}]$. For any
other $r^{\prime}_{1}\in[R_{1}]/\\{r_{1}\\}$, by Jensen’s inequality we have

$\frac{1}{W}[\mathcal{N}\times_{1}\bm{D}^{(1)^{T}}\times_{2}\cdots\times_{K}\bm{D}^{(K)^{T}}]_{r^{\prime}_{1},\ldots,r_{K}}-\mu^{2}_{r^{\prime}_{1},\ldots,r_{K}}\geq
0.$ (37)

Combining the sum of (A.4.2) and (37) over $(r_{2},\ldots,r_{K})$ gives

$\displaystyle
G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})-\sum_{r_{1},\ldots,r_{K}}p^{(1)}_{r_{1}}\cdots
p^{(K)}_{r_{K}}c^{2}_{r_{1},\ldots,r_{K}}$ $\displaystyle\leq$
$\displaystyle-\varepsilon\displaystyle\sum_{r_{2},\ldots,r_{K}}\frac{D_{a_{2}r_{2}}^{(2)}\cdots
D_{a_{K}r_{K}}^{(K)}}{4}(c_{a_{1},a_{2},\ldots,a_{K}}-c_{a_{1}^{\prime},a_{2},\ldots,a_{K}})^{2}$
$\displaystyle\leq$
$\displaystyle-\frac{1}{4}\varepsilon\tau^{K-1}\delta_{\min},$

where the last line uses the fact that
$\displaystyle\sum_{r_{k}}D_{a_{k}r_{k}}^{(k)}=p_{a_{k}}^{(k)}\geq\tau$. ∎

####  A.4.3 Proof of Theorem 2

###### Proof.

The notations we use here are inherited from Lemma 2. By Lemma 2, we obtain
that

$\displaystyle\mathbb{P}\left(\operatorname*{MCR}(\hat{\bm{M}}_{k},\bm{M}_{k,\text{true}})\geq\varepsilon\right)$
$\displaystyle\leq$
$\displaystyle\mathbb{P}\left(G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})-G(\bm{M}_{1,\text{true}},\ldots,\bm{M}_{K,\text{true}})\leq-\frac{1}{4}\varepsilon\tau^{K-1}\delta_{\min}\right).$
(38)

Define
$r=\displaystyle\sup_{\mathcal{J}_{\tau}}|F(\bm{M}_{1},\ldots,\bm{M}_{K})-G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})|$
as the stochastic deviation caused by the label assignment. When the event
$G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})-G(\bm{M}_{1,\text{true}},\ldots,\bm{M}_{K,\text{true}})\leq-\frac{1}{4}\varepsilon\tau^{K-1}\delta_{\min}$
holds, by triangle inequality, we have

$F(\hat{\bm{M}}_{1},\ldots,\hat{\bm{M}}_{K})-F(\bm{M}_{1,\text{true}},\ldots,\bm{M}_{K,\text{true}})\leq
2r-\frac{1}{4}\varepsilon\tau^{K-1}\delta_{\min}.$ (39)

Plugging the event (39) back into inequality (A.4.3), we obtain

$\displaystyle\mathbb{P}\left(\operatorname*{MCR}(\hat{\bm{M}}_{k},\bm{M}_{k,\text{true}})\geq\varepsilon\right)$
$\displaystyle\leq$
$\displaystyle\mathbb{P}\left(F(\hat{\bm{M}}_{1},\ldots,\hat{\bm{M}}_{K})-F(\bm{M}_{1,\text{true}},\ldots,\bm{M}_{K,\text{true}})\leq
2r-\frac{1}{4}\varepsilon\tau^{K-1}\delta_{\min}\right)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(r\geq\frac{\varepsilon\tau^{K-1}\delta_{\min}}{8}\right),$
(40)

where the last line uses the fact that the ${\hat{\bm{M}}_{k}}$ is the global
optimizer of $F(\cdot)$; i.e.
$F(\hat{\bm{M}}_{1},\ldots,\hat{\bm{M}}_{K})=\arg\max
F(\bm{M}_{1},\ldots,\bm{M}_{K})\geq
F(\bm{M}_{1,\text{true}},\ldots,\bm{M}_{K,\text{true}})$.

Now we aim to find the probability (A.4.3) with respect to
$r=\displaystyle\sup_{\mathcal{J}_{\tau}}|F(\bm{M}_{1},\ldots,\bm{M}_{K})-G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})|$.
Note that $r$ involves the quadratic objective $f(x)=x^{2}$. The quadratic
function $f(x)$ is locally lipschitz continuous with lipschitz constant
$b=\sup_{x}|f^{\prime}(x)|$, where $x$ is in the closure of the convex hull of
the entries of $\mathcal{C}$. Note that $b\leq
2\lVert\mathcal{C}\rVert_{\max}$. Therefore, for any partitions
$\\{\bm{M}_{k}\\}$ (which are not necessarily equal to
$\\{\hat{\bm{M}}_{k}\\}$ or $\\{\bm{M}_{k,\text{true}}\\}$):

$\displaystyle\left|F(\bm{M}_{1},\ldots,\bm{M}_{K})-G(\bm{D}^{(1)},\ldots,\bm{D}^{(K)})\right|$
$\displaystyle\leq$
$\displaystyle\sum_{r_{1},\ldots,r_{K}}p_{r_{1}}^{(1)}p_{r_{2}}^{(2)}\cdots
p_{r_{K}}^{(K)}\left|f(\hat{c}_{r_{1},\ldots,r_{K}})-f(\mu_{r_{1},\ldots,r_{K}})\right|$
$\displaystyle\leq$ $\displaystyle\
2\lVert\mathcal{C}\rVert_{\max}\lVert\mathcal{R}(\bm{M}_{1},\ldots,\bm{M}_{K})\rVert_{\max},$
(41)

where

$\hat{c}_{r_{1},\ldots,r_{K}}={1\over\prod_{k}p^{(k)}_{r_{k}}}(\mathcal{Y}\times_{1}\bm{M}^{T}_{1}\times_{2}\cdots\times_{K}\bm{M}^{T}_{K})_{r_{1},\ldots,r_{K}},$

and

$\mu_{r_{1},\ldots,r_{K}}={1\over\prod_{k}p^{(k)}_{r_{k}}}\left[\mathcal{C}\times_{1}\bm{D}^{(1)^{T}}\times_{2}\cdots\times_{K}\bm{D}^{(K)^{T}}\right]_{r_{1},\ldots,r_{K}}.$

are, respectively, sample average and expected sample average, conditional on
the partitions ${\bm{M}_{k}}$, and $\mathcal{R}(\bm{M}_{1},\ldots,\bm{M}_{K})$
is the residual tensor defined in (28).

Combining (A.4.3), (A.4.3) and Hoeffding’s inequality, we have

$\displaystyle\mathbb{P}\left(\operatorname*{MCR}(\hat{\bm{M}}_{k},\bm{M}_{k,\text{true}})\geq\varepsilon\right)$
$\displaystyle\leq\mathbb{P}\left(\sup_{\mathcal{J}_{\tau}}\lVert\mathcal{R}(\bm{M}_{1},\ldots,\bm{M}_{K})\rVert_{\max}\geq{\varepsilon\tau^{K-1}\delta_{\min}\over
16\lVert\mathcal{C}\rVert_{\max}}\right)$
$\displaystyle\leq\mathbb{P}\left(\sup_{I\in\mathcal{I}}{\left|\sum_{(i_{1},\ldots,i_{K})\in
I}\left(Y_{i_{1},\ldots,i_{K}}-\mathbb{E}(Y_{i_{1},\ldots,i_{K}})\right)\right|\over|I|}\geq{\varepsilon\tau^{K-1}\delta_{\min}\over
16\lVert\mathcal{C}\rVert_{\max}}\right)$ $\displaystyle\leq
2^{1+\sum_{k}d_{k}}\text{exp}\left(-\frac{\varepsilon^{2}\tau^{2(K-1)}\delta_{\min}^{2}L}{512\sigma^{2}\lVert\mathcal{C}\rVert_{\max}^{2}}\right),$
(42)

where the last line uses the sub-Gaussianness of the entries in the residual
tensor (conditional on $\\{\bm{M}_{k}\\}$), and
$L=\inf\\{|I|:I\subset\mathcal{I}\\}\geq\tau^{K}\prod_{k=1}^{K}d_{k}$ is
introduced in Section A.4.1. Defining $C=\frac{1}{512}$ in (A.4.3) yields the
desired conclusion. ∎

###  A.5 Sparse estimator

######  Lemma 3.

Consider the regularized least-square estimation,

$\hat{\Theta}^{\text{sparse}}=\operatorname*{arg\,min}_{\Theta\in\mathcal{P}}\left\\{\lVert\mathcal{Y}-\Theta\rVert_{F}^{2}+\lambda\lVert\mathcal{C}\rVert_{\rho}\right\\},$
(43)

where $\mathcal{C}=\llbracket
c_{r_{1},\ldots,r_{K}}\rrbracket\in\mathbb{R}^{R_{1}\times\cdots\times R_{K}}$
is the block-mean tensor, $\lVert\mathcal{C}\rVert_{\rho}$ is the penalty
function with $\rho$ being an index for the tensor norm, and $\lambda$ is the
penalty tuning parameter. We have

$\hat{c}^{\text{sparse}}_{r_{1},\ldots,r_{K}}=\begin{cases}\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}}\mathds{1}\left\\{|\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}}|\geq{\sqrt{\lambda\over
n_{r_{1},\ldots,r_{K}}}}\right\\}&\text{if}\ \rho=0,\\\
\text{sign}(\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}})\left(|\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}}|-{\lambda\over
2n_{r_{1},\ldots,r_{K}}}\right)_{+}&\text{if}\ \rho=1,\end{cases}$ (44)

where $a_{+}=\max(a,0)$ and $\hat{c}^{\text{ols}}_{r_{1},\ldots,r_{K}}$
denotes the ordinary least-square estimate as in Algorithm 1.

###### Proof.

We formulate the estimation of $\mathcal{C}$ as a regularized least-square
regression. Note that $\Theta\in\mathcal{P}$ implies that

$\Theta=\mathcal{C}\times_{1}\bm{M}_{1}\times\cdots\times_{K}\bm{M}_{K}.$ (45)

Define $\bm{X}=\bm{M}_{1}\otimes\ldots\otimes\bm{M}_{K}\in\mathbb{R}^{d\times
R}$, where $d=\prod_{k}d_{k}$ and $R=\prod_{k}R_{k}$, and
${\bm{\beta}}=\text{vec}(\mathcal{C})\in\mathbb{R}^{R}$. Here $\bm{X}$ is a
membership matrix that indicates the block allocation among tensor entries.
Specifically, $\bm{X}$ consists of orthogonal columns with
$\bm{X}^{T}\bm{X}=\text{diag}(n_{1},\ldots,n_{R})$, where $n_{r}$ is the
number of entries in the tensor block that corresponds to the $r$-th column of
$\bm{X}$.

For a given set of $\bm{M}_{k}^{\prime}s$, the optimization (45) with respect
to $\mathcal{C}$ is equivalent to a regularized linear regression with
$\bm{Y}=\text{vec}(\mathcal{Y})$ as the response and $\bm{X}$ as the design
matrix:

$L({\bm{\beta}})=\lVert\bm{Y}-\bm{X}\bm{\beta}\rVert_{2}^{2}+\lambda\lVert\bm{\beta}\rVert_{\rho}.$
(46)

When $\lambda=0$ (no penalty), the minimizer is
$\hat{\bm{\beta}}^{\text{ols}}=(\hat{\beta}^{\text{ols}}_{1},\ldots,\hat{\beta}^{\text{ols}}_{R})=(\bm{X}^{T}\bm{X})^{-1}\bm{X}^{T}\bm{Y}$,
where $\hat{\beta}^{\text{ols}}_{r}={1\over
n_{r}}\bm{y}_{\bm{r}}\mathbf{1}^{T}_{n_{r}}$ for all $r\in[R]$.

Case 1:  $\rho=0$.

Note that $\bm{X}$ induces a partition of indices $[d]$ into $R$ blocks. With
a little abuse of notation, we use $\mathcal{R}=\\{i\in[d]:\bm{X}(i)=r\\}$ to
denote the collection of tensor indices that belong to the $r$th block, and
use $\bm{Y}_{\mathcal{R}}\in\mathbb{R}^{n_{r}}$ to denote the corresponding
tensor entries. By the orthogonality of $\bm{X}$, we have

$\begin{split}L(\bm{\beta})&=\sum_{r=1}^{R}\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda\sum_{r=1}^{R}\mathds{1}{\\{\beta_{r}\neq
0\\}}\\\
&=\sum_{r=1}^{R}\mathop{\mathchoice{\underbrace{\displaystyle\left(\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda\mathds{1}{\\{\beta_{r}\neq
0\\}}\right)}}{\underbrace{\textstyle\left(\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda\mathds{1}{\\{\beta_{r}\neq
0\\}}\right)}}{\underbrace{\scriptstyle\left(\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda\mathds{1}{\\{\beta_{r}\neq
0\\}}\right)}}{\underbrace{\scriptscriptstyle\left(\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda\mathds{1}{\\{\beta_{r}\neq
0\\}}\right)}}}\limits_{:=L_{r}(\beta_{r})}\end{split}$ (47)

The optimization can be separated into each of $\beta_{r}$’s. For any
$r\in[R]$, the sub-optimization $\min_{\beta_{r}}L_{r}(\beta_{r})$ has a
closed-form solution

$\min_{\beta_{r}}L_{r}(\beta_{r})=\begin{cases}\bm{Y}^{T}_{\mathcal{R}}\bm{Y}_{\mathcal{R}}-n_{r}\left(\hat{\beta}_{r}^{\text{ols}}\right)^{2}+\lambda&\text{if
$\hat{\beta}^{\text{ols}}_{r}\neq 0$,}\\\
\bm{Y}^{T}_{\mathcal{R}}\bm{Y}_{\mathcal{R}}&\text{if
$\hat{\beta}^{\text{ols}}_{r}=0$,}\end{cases}$

with

$\arg\min_{\beta_{r}}L_{r}(\beta_{r})=\begin{cases}0\quad&\text{if
$n_{r}\left(\hat{\beta}_{r}^{\text{ols}}\right)^{2}\leq\lambda$},\\\
\hat{\beta}_{r}^{\text{ols}}&\text{otherwise}.\end{cases}$ (48)

Solution (48) can be simplified as
$\hat{\beta}^{\text{sparse}}_{r}=\hat{\beta}^{\text{ols}}_{r}\mathds{1}\\{|\hat{\beta}^{\text{ols}}_{r}|\leq\sqrt{\lambda\over
n_{r}}\\}$. The proof is complete by noting that
$\hat{c}^{\text{sparse}}_{r_{1},\ldots,r_{R}}=\hat{\beta}^{\text{sparse}}_{r}$
and $n_{r_{1},\ldots,r_{K}}=n_{r}$ for all
$(r_{1},\ldots,r_{K})\in[R_{1}]\times\cdots\times[R_{K}]$.

Case 2:  $\rho=1$.

Similar as in Case 1, we write the optimization (46) as

$L({\bm{\beta}})=\sum_{r=1}^{R}\mathop{\mathchoice{\underbrace{\displaystyle\left(\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda|\beta_{r}|\right)}}{\underbrace{\textstyle\left(\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda|\beta_{r}|\right)}}{\underbrace{\scriptstyle\left(\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda|\beta_{r}|\right)}}{\underbrace{\scriptscriptstyle\left(\lVert\bm{Y}_{\mathcal{R}}-\beta_{r}\mathbf{1}_{n_{r}}\rVert_{2}^{2}+\lambda|\beta_{r}|\right)}}}\limits_{:=L_{r}(\beta_{r})},$

where, with a little abuse of notation, we still use $L_{r}(\beta_{r})$ to
denote the sub-optimization. To solve
$\operatorname*{arg\,min}_{\beta_{r}}L_{r}(\beta_{r})$, we use the properties
of subderivative. Taking the subderivative with respect to $\beta_{r}$, we
obtain

$\frac{\partial
L_{r}(\beta_{r})}{\partial\beta_{r}}=\begin{cases}2n_{r}\beta_{r}-2n_{r}\hat{\beta}^{\text{ols}}_{r}+\lambda&\mbox{if
$\beta_{r}>0$,}\\\ [2n_{r}\beta_{r}-2\hat{\beta}^{\text{ols}}_{r}-\lambda,\
2n_{r}\beta_{r}-\hat{\beta}^{\text{ols}}+\lambda]&\mbox{if $\beta_{r}=0$,}\\\
2n_{r}\beta_{r}-2n_{r}\hat{\beta}^{\text{ols}}_{r}+\lambda&\mbox{if
$\beta_{r}<0$.}\\\ \end{cases}$ (49)

Because $\hat{\beta}^{\text{sparse}}_{r}$ minimizes $L_{r}(\beta_{r})$ if and
only if $0\in\frac{\partial L_{r}(\beta_{r})}{\partial\beta_{j}}$, we have:

$\hat{\beta}^{\text{sparse}}_{r}=\begin{cases}\hat{\beta}^{\text{ols}}_{r}+{\lambda\over
2n_{r}}&\mbox{if $\hat{\beta}^{\text{ols}}_{r}<-{\lambda\over 2n_{r}}$,}\\\
0&\mbox{if $\hat{\beta}^{\text{ols}}_{r}\in[-{\lambda\over
2n_{r}},{\lambda\over 2n_{r}}]$,}\\\
\hat{\beta}^{\text{ols}}_{r}-{\lambda\over 2n_{r}}&\mbox{if
$\hat{\beta}^{\text{ols}}_{r}>{\lambda\over 2n_{r}}$.}\\\ \end{cases}$ (50)

The solution (50) can be simplified as

$\hat{\beta}^{\text{sparse}}_{r}=\text{sign}(\hat{\beta}^{\text{ols}}_{r})\left(|\hat{\beta}^{\text{ols}}_{r}|-{\lambda\over
2n_{r}}\right)_{+},\quad\text{for all $r\in[R]$}.$

∎

##  Appendix B Supplementary Figures and Tables

Supplementary Figure S1: (a) estimation error and (b) sparse error rate against noise for sparse tensors of dimension $(40,40,40)$ when $p=0.8$.  Dimensions | True clustering sizes | Noise | Estimated clustering sizes
---|---|---|---
$(d_{1},d_{2},d_{3})$ | $(R_{1},R_{2},R_{3})$ | $(\sigma)$ | $(\hat{R}_{1},\hat{R}_{2},\hat{R}_{3})$
$(40,40,40)$ | $(4,4,4)$ | 4 | $({\bf 4},\ {\bf 4},\ {\bf 4})\pm(0,\ 0,\ 0)$
$(40,40,40)$ | $(4,4,4)$ | 8 | $({\bf 3.94},\ {\bf 3.96},\ {\bf 3.96})\pm(0.03,\ 0.03,\ 0.03)$
$(40,40,40)$ | $(4,4,4)$ | 12 | $(3.08,\ 3.12,\ 3.12)\pm(0.10,0.10,0.10)$
$(40,40,80)$ | $(4,4,4)$ | 4 | $({\bf 4},\ {\bf 4},\ {\bf 4})\pm(0,\ 0,\ 0)$
$(40,40,80)$ | $(4,4,4)$ | 8 | $({\bf 4},\ {\bf 4},\ {\bf 4})\pm(0,\ 0,\ 0)$
$(40,40,80)$ | $(4,4,4)$ | 12 | $({\bf 3.96},\ {\bf 3.96},\ 3.92)\pm(0.04,0.04,0.04)$
$(40,40,40)$ | $(2,3,4)$ | 4 | $({\bf 2},\ {\bf 3},\ {\bf 4})\pm(0,\ 0,\ 0)$
$(40,40,40)$ | $(2,3,4)$ | 8 | $({\bf 2},{\bf\ 3},\ {\bf 3.96})\pm(0,\ 0,\ 0.03)$
$(40,40,40)$ | $(2,3,4)$ | 12 | $({\bf 2},\ {\bf 2.96},\ 3.60)\pm(0,\ 0.05,\ 0.09)$
Supplementary Table S1: The simulation results for estimating
$\bm{R}=(R_{1},R_{2},R_{3})$. Bold number indicates no significant difference
between the estimate and the ground truth, based on a $z$-test with a level
$0.05$.

Tissues | Over-expressed genes | Block-means | Under-expressed genes | Block-means
---|---|---|---|---
Cluster 1 | GFAP, MBP | 10.88 | GPR6 , DLX5 , DLX6 , NKX2-1 | -8.40
Cluster 2 | GFAP, MBP | 5.98 | CDH9, RXFP1, CRH, ARX, CARTPT, DLX1,FEZF2 | -9.49
Cluster 3 | GFAP, MBP | 8.34 | AVPR1A, CCKAR, CHRNB4, CYP19A1, HOXA4 , LBX1, SLC6A3 | -8.45
TBR1, SLC17A6, SLC30A3 | -8.17
Cluster 4 | GFAP, MBP | 8.83 | AVPR1A, CCKAR, CHRNB4, CYP19A1, HOXA4 , LBX1, SLC6A3 | -8.40
DAO EN2 EOMES | -6.57

Supplementary Table S2: Top expression blocks from the multi-tissue gene
expression analysis. The tissue clusters are described in Supplementary
Section D.

Countries | Countries | Relation types
---|---|---
Cluster 1 | Clusters 4 and 5 | reltreaties, booktranslations, relbooktranslations, relexports, exports3
Clusters 1 and 4 | Cluster 5 | relintergovorgs, relngo, intergovorgs3, ngoorgs3
Cluster 3 | Clusters 1, 4, and 5 | commonbloc0, blockpositionindex
Clusters 1 and 3 | Clusters 4 and 5 | timesinceally, independence
Cluster 1 | Cluster 3
Cluster 4 | Cluster 5
Cluster 4 | Cluster 5 | treaties, conferences, weightedunvote, unweightedunvote, intergovorgs, ngo,
officialvisits, exportbooks, relexportbooks, tourism,
reltourism, tourism3, exports, militaryalliance, commonbloc2

Supplementary Table S3: Top blocks from the _Nations_ data analysis. The
countries clusters are described in Supplementary Section D.

##  Appendix C Time complexity

The total cost of our Algorithm 1 is $\mathcal{O}(d)$ per iteration, where
$d=\prod_{k}d_{k}$ denotes the total number of tensor entries. The per-
iteration computational cost scales linearly with the sample size, and this
complexity is comparable to the classical tensor methods such as CP and Tucker
decomposition. More specifically, each iteration of Algorithm 1 consists of
updating the core tensor $\mathcal{C}$ and $K$ membership matrices
$\bm{M}_{k}$’s. The update of $\mathcal{C}$ requires $\mathcal{O}(d)$
operations and the update of $\bm{M}_{k}$ requires $\mathcal{O}(R_{k}{d\over
d_{k}})$ operations. Therefore the total cost is
$\mathcal{O}(d+d\sum_{k}{R_{k}\over d_{k}})$.

##  Appendix D Additional information for real data analysis

Multi-tissue gene expression. The gene expression data we analyzed is part of
the GTEx v6 datasets (https://www.gtexportal.org/home/datasets). We cleaned
and preprocessed the data following the steps in [4]. We focused on the 13
brain tissues, 193 individuals, and 362 annotated genes provided by Atlax of
the Developing Human Brain (http://www.brainspan.org/ish). After applying the
$\ell$-0 penalized TBM to the mean-centered data tensor, we identified the
following four clusters of tissues:

  1. -

Cluster 1: Substantia nigra, Spinal cord (cervical c-1)

  2. -

Cluster 2: Cerebellum, Cerebellar Hemisphere

  3. -

Cluster 3: Caudate (basal ganglia), Nucleus accumbens (basal ganglia), Putamen
(basal ganglia)

  4. -

Cluster 4: Cortex, Hippocampus, Anterior cingulate cortex (BA24), Frontal
Cortex (BA9), Hypothalamus, Amygdala

We found that most tissue clusters are spatially restricted to specific brain
regions, such as the two cerebellum tissues (cluster 2), three basal ganglia
tissues (cluster 3), and the cortex tissues (cluster 4). Supplementary Table
S2 reports the associated gene cluster for each tissue cluster. Because our
method attaches importance to blocks by the absolute mean estimates, our
method is able to detect both over- and under-expression patterns. Blocks with
highly positive means correspond to over-expressed genes, whereas blocks with
highly negative means correspond to under-expressed genes.

Nations dataset. This is a $14\times 14\times 56$ binary tensor consisting of
$56$ political relations of $14$ countries between 1950 and 1965 [2]. The
tensor entry indicates the presence or absence of a political action, such as
“treaties”, “sends tourists to”, between the nations. We applied the $\ell$-0
penalized TBM to the binary-valued data tensor, and we identified the
following five clusters of countries:

  1. -

Cluster 1: Brazil, Egypt, India, Israel, Netherlands

  2. -

Cluster 2: Burma, Indonesia, Jordan

  3. -

Cluster 3: China, Cuba, Poland, USSA

  4. -

Cluster 4: USA

  5. -

Cluster 5: UK

Supplementary Table S3 reports the cluster constitutions for top blocks.
Because the tensor entries take value on either 0 or 1, the top blocks mostly
have mean one.

## References

  * [1] Hua Zhou, Lexin Li, and Hongtu Zhu.  Tensor regression with applications in neuroimaging data analysis.  Journal of the American Statistical Association, 108(502):540–552, 2013.
  * [2] Maximilian Nickel, Volker Tresp, and Hans-Peter Kriegel.  A three-way model for collective learning on multi-relational data.  In International Conference on Machine Learning, volume 11, pages 809–816, 2011.
  * [3] Yichuan Tang, Ruslan Salakhutdinov, and Geoffrey Hinton.  Tensor analyzers.  In International Conference on Machine Learning, pages 163–171, 2013.
  * [4] Miaoyan Wang, Jonathan Fischer, and Yun S Song.  Three-way clustering of multi-tissue multi-individual gene expression data using constrained tensor decomposition.  Annals of Applied Statistics, in press, 2019.
  * [5] Victoria Hore, Ana Viñuela, Alfonso Buil, Julian Knight, Mark I McCarthy, Kerrin Small, and Jonathan Marchini.  Tensor decomposition for multiple-tissue gene expression experiments.  Nature genetics, 48(9):1094, 2016.
  * [6] Frank L Hitchcock.  The expression of a tensor or a polyadic as a sum of products.  Journal of Mathematics and Physics, 6(1-4):164–189, 1927.
  * [7] Ledyard R Tucker.  Some mathematical notes on three-mode factor analysis.  Psychometrika, 31(3):279–311, 1966.
  * [8] Kean Ming Tan and Daniela M Witten.  Sparse biclustering of transposable data.  Journal of Computational and Graphical Statistics, 23(4):985–1008, 2014.
  * [9] Tamara G Kolda and Jimeng Sun.  Scalable tensor decompositions for multi-aspect data mining.  In 2008 Eighth IEEE international conference on data mining, pages 363–372. IEEE, 2008.
  * [10] Chang-Dong Wang, Jian-Huang Lai, and S Yu Philip.  Multi-view clustering based on belief propagation.  IEEE Transactions on Knowledge and Data Engineering, 28(4):1007–1021, 2015.
  * [11] Miaoyan Wang and Lexin Li.  Learning from binary multiway data: Probabilistic tensor decomposition and its statistical optimality.  arXiv preprint arXiv:1811.05076, 2018.
  * [12] Stefanie Jegelka, Suvrit Sra, and Arindam Banerjee.  Approximation algorithms for tensor clustering.  In International Conference on Algorithmic Learning Theory, pages 368–383. Springer, 2009.
  * [13] Eric C Chi, Brian R Gaines, Will Wei Sun, Hua Zhou, and Jian Yang.  Provable convex co-clustering of tensors.  arXiv preprint arXiv:1803.06518, 2018.
  * [14] Tamara G Kolda and Brett W Bader.  Tensor decompositions and applications.  SIAM review, 51(3):455–500, 2009.
  * [15] Sara C Madeira and Arlindo L Oliveira.  Biclustering algorithms for biological data analysis: a survey.  IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), 1(1):24–45, 2004.
  * [16] Chao Gao and Zongming Ma.  Minimax rates in network analysis: Graphon estimation, community detection and hypothesis testing.  arXiv preprint arXiv:1811.06055, 2018.
  * [17] Peter J Bickel and Aiyou Chen.  A nonparametric view of network models and newman–girvan and other modularities.  Proceedings of the National Academy of Sciences, 106(50):21068–21073, 2009.
  * [18] Kehui Chen Jing Lei and Brian Lynch.  Consistent community detection in multi-layer network data.  Biometrika, to appear, 2019.
  * [19] Subhadeep Paul, Yuguo Chen, et al.  Consistent community detection in multi-relational data through restricted multi-layer stochastic blockmodel.  Electronic Journal of Statistics, 10(2):3807–3870, 2016.
  * [20] Robin A Darton.  Rotation in factor analysis.  Journal of the Royal Statistical Society: Series D (The Statistician), 29(3):167–194, 1980.
  * [21] Hervé Abdi.  Factor rotations in factor analyses.  Encyclopedia for Research Methods for the Social Sciences, Sage: Thousand Oaks, pages 792–795, 2003.
  * [22] Anru Zhang and Dong Xia.  Tensor SVD: Statistical and computational limits.  IEEE Transactions on Information Theory, 2018.
  * [23] Chao Gao, Yu Lu, Zongming Ma, and Harrison H Zhou.  Optimal estimation and completion of matrices with biclustering structures.  The Journal of Machine Learning Research, 17(1):5602–5630, 2016\.
  * [24] Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M Kakade, and Matus Telgarsky.  Tensor decompositions for learning latent variable models.  The Journal of Machine Learning Research, 15(1):2773–2832, 2014\.
  * [25] Miaoyan Wang and Yun Song.  Tensor decompositions via two-mode higher-order SVD (HOSVD).  In Artificial Intelligence and Statistics, pages 614–622, 2017\.
  * [26] Daniel Aloise, Amit Deshpande, Pierre Hansen, and Preyas Popat.  NP-hardness of Euclidean sum-of-squares clustering.  Machine learning, 75(2):245–248, 2009.
  * [27] Nuala A O’Leary, Mathew W Wright, J Rodney Brister, Stacy Ciufo, Diana Haddad, Rich McVeigh, Bhanu Rajput, Barbara Robbertse, Brian Smith-White, Danso Ako-Adjei, et al.  Reference sequence (refseq) database at ncbi: current status, taxonomic expansion, and functional annotation.  Nucleic acids research, 44(D1):D733–D745, 2015.
  * [28] Phillippe Rigollet and Jan-Christian Hütter.  High dimensional statistics.  Lecture notes for course 18S997, 2015.