11footnotetext: Department of Statistics, University of Wisconsin-
Madison22footnotetext: Department of Statistics, Columbia
University33footnotetext: Department of Biostatistics & Bioinformatics, Duke
University

# A Sharp Blockwise Tensor Perturbation Bound for Orthogonal Iteration

Yuetian Luo1, Garvesh Raskutti1, Ming Yuan2, and Anru R. Zhang1,3

###### Abstract

In this paper, we develop novel perturbation bounds for the high-order
orthogonal iteration (HOOI) [DLDMV00b]. Under mild regularity conditions, we
establish blockwise tensor perturbation bounds for HOOI with guarantees for
both tensor reconstruction in Hilbert-Schmidt norm
$\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm HS}}$ and mode-$k$
singular subspace estimation in Schatten-$q$ norm
$\|\sin\Theta(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k})\|_{q}$ for any
$q\geq 1$. We show the upper bounds of mode-$k$ singular subspace estimation
are unilateral and converge linearly to a quantity characterized by blockwise
errors of the perturbation and signal strength. For the tensor reconstruction
error bound, we express the bound through a simple quantity $\xi$, which
depends only on perturbation and the multilinear rank of the underlying
signal. Rate matching deterministic lower bound for tensor reconstruction,
which demonstrates the optimality of HOOI, is also provided. Furthermore, we
prove that one-step HOOI (i.e., HOOI with only a single iteration) is also
optimal in terms of tensor reconstruction and can be used to lower the
computational cost. The perturbation results are also extended to the case
that only partial modes of ${\mathbfcal{T}}$ have low-rank structure. We
support our theoretical results by extensive numerical studies. Finally, we
apply the novel perturbation bounds of HOOI on two applications, tensor
denoising and tensor co-clustering, from machine learning and statistics,
which demonstrates the superiority of the new perturbation results.

###### Abstract

In this supplement, we provide an additional corollary of our main results and
all proofs in this paper.

##  1 Introduction

The past decades have seen a large body of work on tensors or multiway arrays
[KB09, SDLF+17, CMDL+15, Kro08]. Tensors arise in numerous applications
involving multiway data (e.g., brain imaging [ZLZ13], hyperspectral imaging
[LL10], recommender system design [BQS18]). In addition, various methods have
been developed and applied to fundamental tensor problems such as tensor
completion [YZ16, XYZ17, YZ17, Zha19], tensor regression [ZLRY20, ZLZ13,
RYC19, CRY16], tensor PCA/SVD [ZX18, RM14, LYZ17], generalized tensor learning
[HWZ20]. In many other problems where the observations are not necessarily
tensors, tensor parameters arise, such as topic and latent variable models
[AGH+14], additive index models [BFY18], high-order interaction pursuit
[HZC18]. We refer readers to recent survey papers [KB09, SDLF+17, CMDL+15].

Among these methods, tensor decomposition is one of the most important and a
flurry of research have been devoted to it. Tensor decomposition plays a
similar role to singular value decomposition (SVD) or eigendecomposition for
matrices which is of fundamental importance throughout a wide range of fields
including computer science, applied mathematics, machine leaning, statistics,
signal processing, etc. For matrices, truncated SVD achieves the best low-rank
approximation in terms of any unitarily invariant norm by the well-known
Eckart-Young-Mirsky theorem [EY36, Mir60, GHS87] and more importantly it is
computationally efficient. Despite the well-established theory for low-rank
decomposition of matrices, tensors present unique challenges. First there are
several notions of low-rankness in tensors, moreover it has been shown that
computing various best low-rank approximations of a given tensor is NP hard in
general [HL13].

Fortunately, many computationally efficient algorithms have been proposed to
approximate the best low-rank tensor decomposition [KDL80, DLDMV00b, ES09,
IAVHDL11, IDLAVH09, SL10, DLDMV00a, VVM12]. One popular choice is the high-
order orthogonal iteration (HOOI) proposed in [DLDMV00b]. HOOI is based on
alternating least-squares. It can be seen as a “spectral” algorithm for
tensors, generalizations of the 2D-PCA [SS07], and the power iteration
refinement of HOSVD [DLDMV00a] and sequential HOSVD [VVM12]. Convergence
properties of HOOI have been studied in [ZG01, WC14, Usc15, Xu14, Xu18, ZX18].

In addition to computing low-rank approximations of matrices and tensors,
there is the more nuanced question of computing low-rank approximations under
noise perturbation and determining how the perturbation impacts the quality of
the decomposition. For matrices, perturbation theory is well studied and a
number of results has been established [Bha13, Ste90]. However, perturbation
theory for tensors is still in its infancy. It is difficult to extend results
for matrices to tensors due to the complexity of tensor algebra and the fact
many well established theories and concepts in matrices such as SVD or
eigendecomposition do not exist or are not easy to use for tensors. There are
several attempts in the literature and most of them require the noise tensor
to be random [RM14, ZX18, HSS15, LYZ17, ZH18, ZT15, BCS+16, XZ19]. In this
paper, we move one step forward in this direction and provide the first
general perturbation bounds of HOOI for tensors under the signal-plus-noise
model

$\widetilde{{\mathbfcal{T}}}={\mathbfcal{T}}+{\mathbfcal{Z}}$ (1)

without putting any structural assumption on the perturbation noise
${\mathbfcal{Z}}$. Like the classical Wedin’s perturbation theory for
matrices, we provide perturbation guarantees of estimated mode-$k$ singular
subspace in tensors. In addition, we also provide the perturbation bounds for
tensor reconstruction. By providing the deterministic rate matching lower
bound, we can further show HOOI with good initialization is optimal for tensor
reconstruction.

###  1.1 Problem Statement

Formally, this paper considers the tensor perturbation model (1), where
${\mathbfcal{T}}$ is the low-rank order-$d$ signal tensor and
${\mathbfcal{Z}}$ is the perturbation tensor with the same dimension as
${\mathbfcal{T}}$. Two popular choices of low-rankness in tensors are
canonical polyadic (CP) low-rank and multilinear/Tucker low-rank and each of
them has their respective benefits (see the surveys [KB09, SDLF+17, CMDL+15,
GKT13]). The CP low-rank decomposition which approximates the original tensor
by a sum of rank-1 outer products gives a compact and unique (under certain
conditions) representation and multilinear/Tucker low-rank decomposition
generally finds a better fit for the data by estimating the subspace of each
mode. Since any CP low-rank tensor can be written as a multilinear low-rank
tensor with a diagonal core tensor, we focus on the latter (Tucker rank) in
this paper. Specifically, we assume ${\mathbfcal{T}}$ admits the following low
multilinear rank (Tucker) decomposition:

${\mathbfcal{T}}={\mathbfcal{S}}\times_{\Omega_{1}}{\mathbf{U}}_{1}\times\cdots\times_{\Omega_{m}}{\mathbf{U}}_{m}.$
(2)

Here, ${\mathbfcal{S}}$ is an order-$d$ core tensor;
$\\{{\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}}\\}_{i=1}^{m}$ are group-$i$
singular subspaces; $\\{\Omega_{i}\\}_{i=1}^{m}$ are symmetric index groups;
and “$\times_{\Omega_{i}}$” is the tensor matrix product along modes in
$\Omega_{i}$. The formal definition of Tucker decomposition and tensor matrix
product will be given in Section 2.

The symmetric index groups $\\{\Omega_{i}\\}_{i=1}^{m}$ represent the specific
symmetricity structure of ${\mathbfcal{T}}$ and satisfy
$\Omega_{i}\subseteq[d],\Omega_{i}\neq\emptyset$,
$\Omega_{i}\bigcap\Omega_{j}=\emptyset$ for $i\neq j$, and
$\bigcup_{i=1}^{m}\Omega_{i}=[d]$. It means by fixing indices outside the
group and any permutation of indices within the group does not change the
value of tensor corresponding to those indices. For example, if
$\Omega_{1}=\\{1,2,\ldots,k\\}$, then fixing coordinates in
$\\{\Omega_{i}\\}_{i=2}^{m}$ i.e., coordinates $k+1,\ldots,d$ and for any
permutation $\delta$ of $[k]$, we have
${\mathbfcal{T}}_{[i_{\delta(1)},\ldots,i_{\delta(k)},i_{k+1},\ldots,i_{d}]}={\mathbfcal{T}}_{[i_{1},\ldots,i_{d}]}$.
In addition to ${\mathbfcal{T}}$, we also assume ${\mathbfcal{S}}$ and
${\mathbfcal{Z}}$ have the same symmetric structure characterized by
$\\{\Omega_{i}\\}_{i=1}^{m}$. The symmetric index groups have two extreme
cases:

  * •

Asymmetric: $\Omega_{i}=\\{i\\}$ for $i=1,\ldots,d$,

  * •

Supersymmetric: $\Omega_{1}=\\{1,2,\ldots,d\\}$.

To simplify the notation, if mode $j\in\Omega_{i}$, then we denote
$j^{\prime}:=i$ as the group index of mode $j$. For symmetric group $i$, the
dimension of ${\mathbfcal{T}}$ in this group is denoted as $p_{i}$ and its
mode rank is denoted as $r_{i}$. The rigorous definition of mode rank is
postponed to Section 2. Also throughout the paper, for $i\in[m]$, let
$\bar{i}=\min\\{j:j\in\Omega_{i}\\}$ be the smallest index in $\Omega_{i}$,
and denote
$\widecheck{\Omega}_{i}:=\Omega_{i}\setminus\\{\bar{i}\\},\underline{\Omega}_{i}:=\bigcup_{j=1}^{i-1}\Omega_{j},\overline{\Omega}_{i}:=\bigcup_{j=i+1}^{d}\Omega_{j}$
with $\underline{\Omega}_{1}=\overline{\Omega}_{d}=\emptyset$. Finally, we
define the signal strength of ${\mathbfcal{T}}$ be
$\lambda=\min_{i=1,\ldots,m}\sigma_{r_{i}}(\mathcal{M}_{\bar{i}}({\mathbfcal{T}}))$,
the smallest singular values of matricization of ${\mathbfcal{T}}$ in modes
$\\{\bar{i}\\}_{i=1}^{m}$.

In summary, for the dimensions of the perturbation model (1), we have
$\widetilde{{\mathbfcal{T}}},{\mathbfcal{T}},{\mathbfcal{Z}}\in\mathbb{R}^{p_{1^{\prime}}\times\cdots\times
p_{d^{\prime}}}$ with symmetric index groups $\\{\Omega_{i}\\}_{i=1}^{m}$. The
HOOI algorithm we study is provided in Algorithm 1. It is worth noting the
original HOOI algorithm in [DLDMV00b] mainly focuses on asymmetric tensor
decomposition and we generalize it to accommodate arbitrary symmetric
structures of ${\mathbfcal{T}}$ characterized by $\\{\Omega_{i}\\}_{i=1}^{m}$.
In addition, in the literature [DLDMV00b, KB09], HOOI often refers to the
overall procedure including both the initialization
$\widetilde{\mathbf{U}}_{i}^{(0)}={\rm
SVD}_{r_{i}}(\mathcal{M}_{i}(\widetilde{{\mathbfcal{T}}}))$ and the orthogonal
iteration updates as detailed in Algorithm 1. We shall point out that this
paper studies the orthogonal iteration with any initializers satisfying some
mild conditions. Our results accommodate different types of initialization
scheme and can be applied to a wide range of scenarios (see Remark 1). We
conclude the algorithm by remarking that the matricization mode we choose in
group $i$ to perform SVD in (3) does not matter due to symmetry. For
simplicity, we choose the smallest index in group $i$.

Algorithm 1 High-Order Orthogonal Iteration for Tensor Decomposition

Input:
$\widetilde{{\mathbfcal{T}}}\in\mathbb{R}^{p_{1^{\prime}}\times\cdots\times
p_{d^{\prime}}}$, symmetric index groups $\\{\Omega_{i}\\}_{i=1}^{m}$,
initialization $\\{\widetilde{{\mathbf{U}}}_{i}^{(0)}\\}_{i=1}^{m}$ with
$\widetilde{{\mathbf{U}}}_{i}^{(0)}\in\mathbb{O}_{p_{i},r_{i}}$, maximum
number of iterations $t_{\max}$.
Output:
$\\{\widehat{{\mathbf{U}}}_{i}\\}_{i=1}^{m},\widehat{{\mathbfcal{T}}}$.

1:For $t=1,\ldots,t_{\max}$, do

  1. (a)

For $i=1$ to $m$, update

$\begin{split}\widetilde{{\mathbf{U}}}_{i}^{(t+1)}={\rm
SVD}_{r_{i}}&\big{(}\mathcal{M}_{\bar{i}}(\widetilde{{\mathbfcal{T}}}\times_{\Omega_{1}}(\widetilde{{\mathbf{U}}}_{1}^{(t+1)})^{\top}\times\cdots\\\
&\times_{\Omega_{i-1}}(\widetilde{{\mathbf{U}}}_{i-1}^{(t+1)})^{\top}\times_{\widecheck{\Omega}_{i}}(\widetilde{{\mathbf{U}}}_{i}^{(t)})^{\top}\times\cdots\times_{\Omega_{d}}(\widetilde{{\mathbf{U}}}_{d}^{(t)})^{\top})\big{)}.\end{split}$
(3)

2:Let $\widehat{{\mathbf{U}}}_{i}=\widetilde{{\mathbf{U}}}_{i}^{(t_{\max})}$
for $i=1,\ldots,m$ and compute
$\widehat{{\mathbfcal{T}}}=\widetilde{{\mathbfcal{T}}}\times_{\Omega_{1}}P_{\widehat{{\mathbf{U}}}_{1}}\times\cdots\times_{\Omega_{m}}P_{\widehat{{\mathbf{U}}}_{m}}.$

###  1.2 Our Contributions

Under the perturbation model (1) and the HOOI Algorithm 1, we make the
following major contributions to the tensor perturbation theory of low-rank
tensor decomposition based on HOOI.

  * •

We provide the first sharp blockwise perturbation bounds of HOOI for tensors
with guarantees for both the estimated mode-$k$ singular subspace and the
tensor reconstruction induced by low-rank approximation of the noise corrupted
tensor. Specifically, in Theorem 3, we show that under suitable initialization
$\\{\widetilde{{\mathbf{U}}}^{(0)}_{i}\\}_{i=1}^{m}$ in Algorithm 1, the upper
bound of
$\max_{i}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{i},{\mathbf{U}}_{i}\right)\right\|_{q}$
converges linearly with respect to the iteration number to a fixed quantity
characterized by ${\mathbfcal{Z}}$ and the signal strength. We also give the
first high-order unilateral perturbation bounds for mode-$k$ singular
subspace. More discussion about this can be found in Remark 3. In addition, a
practically useful tensor reconstruction error bound in Hilbert-Schmidt norm
is provided. Surprisingly, we found the upper bound of tensor reconstruction
is free of the “condition number” of the signal tensor and can be expressed by
a unified simple quantity $\xi$ characterized only by the noise tensor
${\mathbfcal{Z}}$ and the underlying multilinear rank of ${\mathbfcal{T}}$.
$\xi$ is closely related to the Gaussian width [Gor88], a common measure for
the complexity of a given set, and its formal definition and more explanations
will be given in Sections 2.1 and 6. It is worth mentioning that the
reconstruction error bound of $\widehat{{\mathbfcal{T}}}$ is often
significantly better than the simple estimator $\widetilde{{\mathbfcal{T}}}$
and initialization estimator truncated HOSVD without power iteration. This is
fundamentally different from the perturbation results in [AGH+14]. See further
discussions in Remarks 4 and 7. In addition, we also generalize the main
results to the case that only a subset of modes of ${\mathbfcal{T}}$ have low-
rank structure.

  * •

Furthermore, we provide a deterministic minimax lower bound for the tensor
reconstruction error under perturbation model (1) in Theorem 2. The lower
bound matches the perturbation upper bound in Theorem 3 when the tensor order
$d$ is fixed, which demonstrates the optimality of HOOI for tensor
reconstruction.

  * •

In addition, by combining Theorems 2 and 3, we prove that the tensor
reconstruction error rate of HOOI with only one iteration is also optimal and
further iterations improve the coefficient in front of the error rate $\xi$.
This suggests that in some applications where running full HOOI is too
expensive and prohibitive compared to truncated HOSVD [DLDMV00a] or
sequentially truncated HOSVD [VVM12], we can just run HOOI for one iteration
to obtain an optimal (up to constant) reconstruction. Details are provided in
Remark 5 and numerical comparison is given in Section 7.2.

  * •

Finally, we perform extensive numerical studies to support our perturbation
bounds and do a comparison with other existing low-rank tensor decomposition
algorithms. In addition, we apply the new perturbation bounds of HOOI in two
modern applications, tensor denoising and tensor co-clustering, from machine
learning and statistics. Based on our perturbation results, we can easily
recover the results of tensor denoising in current literature with a much
shorter proof and provide the first guarantee of HOOI on tensor cocluster
recovery that improves the state-of-the-art results.

###  1.3 Related Literature

In this section, we give a brief overview of the literature on Tucker
decomposition of tensors, matrix/tensor perturbation theory related to this
article. Tensor decomposition has become one of the most important topics in
the literature on tensors [KB09, SDLF+17, CMDL+15]. An analogy of matrix SVD
for Tucker decomposition of tensors, which is today commonly called high-order
singular value decomposition (HOSVD), was first proposed in [Tuc66] and then
popularized by [DLDMV00a]. However, unlike the truncated matrix SVD, truncated
HOSVD (T-HOSVD) can provide a reasonable but not necessarily optimal low-
multilinear-rank approximation for a given tensor. It has been shown that
computing the best low multilinear rank approximation is NP hard in general
[HL13]. On the other hand, various computationally efficient algorithms have
been proposed to obtain better approximations than HOSVD [DLDMV00b, ES09,
IAVHDL11, IDLAVH09, SL10]. Among them, the high-order orthogonal iteration
(HOOI) in [DLDMV00b] has become a popular choice in literature. HOOI has been
applied to problems including but not limited to tensor PCA/tensor SVD [RM14,
ZX18], tensor completion [XYZ17], tensor regression [ZLRY20], hypergraph
community recovery [KSX19, JLLX20], independent component analysis [DLV04],
multiple factor analysis [VT03], tensor clustering [LZ20b]. Many variants of
HOOI such as sparse high-order singular value decomposition (STAT-SVD) for
tensors [ZH18], regularized HOOI [KSX19, JLLX20], generalized higher-order
orthogonal iteration (gHOI) [LSF+14] have been proposed. Nowadays, HOOI has
become a prevalent choice to compute the low-rank tensor approximation in many
applications and been coded in common tensor software such as Matlab “Tensor
Toolbox” [BK12], “Tensorlab” [SVBDL14] and R “rTensor” package [LBW18].
Moreover, it has been regarded as the gold standard guideline for comparison
when developing even faster randomized or memory-efficient algorithms for low-
rank tensor approximation [SGL+19, MB18, CC10, Tso10, KS08].

Perturbation theory is a long-existing field in mathematics. In particular,
the perturbation theory on matrices has attracted much attention. The original
work in matrix dated back to Weyl, Davis-Kahan and Wedin [DK70, Wed72, Wey12,
Ste98] (see [Bha13, Ste90] for an overview of classical perturbation results
and historical development) and recently it has been further developed in
[YWS14, CZ18, CTP19b]. In addition, various generalizations and extensions
have been made in different settings including random perturbation [Vu11,
SN13, OVW18, Wan15, BGN11, AFWZ17, KX16, BGN12, CTP19a, MSC20, CCF18],
structured perturbation [FWZ16, Ste06] and many others [Dop00, EBW17]. Also
the perturbation theory for matrices has been widely applied to a number of
applications such as community detection [RCY11, CCT12, CRV15, STFP12,
CTP19b], covariance matrix estimation [FWZ16, CTP19b], matrix denoising
[CZ18], matrix completion [CCZ16], etc.

Studying perturbation theory for low-rank tensor approximation is much harder
than matrices and there are only a few attempts in literature. Perturbation
results for the best low multilinear rank approximation have been developed in
[DL04, ES11]. These results are promising but less practically useful due to
the computational hardness of the best multilinear rank approximation [HL13].
Hence more attention are given to the perturbation results for polynomial-time
low-rank tensor approximation algorithms. A line of work has been done when
the perturbation is random [RM14, ZX18, HSS15, LYZ17, ZH18, ZT15, BCS+16,
XZ19] and various perturbation results for robust tensor decomposition with
sparse noise have been studied in [GQ14, GGH14, SRT15, AJSN16]. Perturbation
bounds for orthogonally decomposable tensors have been studied in [AY20]. In
addition, [AGH+14, AGJ14] provided perturbation guarantees for power iteration
algorithm for symmetric orthogonal and non-orthogonal CP low-rank
decomposition. However, we are not aware of any perturbation result for
polynomial time algorithms under partial symmetric multilinear low-rank
setting. In this paper, we make the first attempt in this direction and
provide the first perturbation bounds of HOOI for tensors in the general
setting. We end this section by remarking that in most situations there is a
trade-off about the quality of low-rank tensor decomposition and computational
cost of the algorithm. For example, computing truncated HOSVD and sequentially
truncated HOSVD (ST-HOSVD) [VVM12] may be much faster than iterative
algorithms such as HOOI, (quasi)Newton-Grassmann method [ES09, SL10],
geometric Newton method [IDLAVH09] and Riemannian trust region scheme
[IAVHDL11] in the large scale settings. On the other hand, these iterative
algorithms achieve higher accuracy. In the perturbation model (1), we show
HOOI could achieve optimal tensor reconstruction error, which is not true for
HOSVD and ST-HOSVD in general. Another contribution of this paper is we show
that one-step HOOI (i.e., $t_{\max}=1$ in Algorithm 1) is also optimal in
terms of tensor reconstruction and could be used as a surrogate of HOOI in
large scale tensor decomposition settings.

###  1.4 Organization of the Paper

The remainder of the article is organized as follows. In Section 2, after a
brief introduction of notation and preliminaries, we define various blockwise
errors of ${\mathbfcal{Z}}$ as the key quantities in our perturbation bounds.
We illustrate our main perturbation theorem in asymmetric order-$3$ case in
Section 3 and at the end of the same section, we provide the deterministic
lower bound for tensor reconstruction. In Section 4, we provide the
perturbation bounds of HOOI applying on a corrupted general partial symmetric
order-$d$ tensor. In Section 5, we discuss the tensor perturbation bounds when
the target tensor has the low-rank structure only along a subset of modes. In
Section 6, we apply our perturbation bounds to two applications, tensor
denoising and tensor co-clustering. In Section 7, we corroborate our
theoretical results by extensive numerical studies. Conclusion and discussions
are provided in Section 8. Due to space constraints, all technical proofs are
postponed to the Appendix.

##  2 Notation and Preliminaries

The following notation will be used throughout this article. For any non-
negative integer $n$, let $[n]=\\{1,\ldots,n\\}$. Lowercase letters (e.g.,
$a,b$), lowercase boldface letters (e.g., ${\mathbf{u}},{\mathbf{v}}$),
uppercase boldface letters (e.g., ${\mathbf{U}},{\mathbf{V}}$), and boldface
calligraphic letters (e.g., ${\mathbfcal{T}},{\mathbfcal{Z}}$) are used to
denote scalars, vectors, matrices, and order-3-or-higher tensors,
respectively. For any two series of numbers, say $\\{a_{n}\\}$ and
$\\{b_{n}\\}$, denote $a\asymp b$ or $a=O(b)$ if $ca_{n}\leq b_{n}\leq Ca_{n}$
or $a_{n}\leq Cb_{n}$ for some uniform constants $c,C>0$. For any matrix
${\mathbf{D}}\in\mathbb{R}^{m\times n}$ with singular value decomposition
$\sum_{i=1}^{m\land
n}\sigma_{i}({\mathbf{D}}){\mathbf{u}}_{i}{\mathbf{v}}_{i}^{\top}$, let
${\mathbf{D}}_{\max(r)}=\sum_{i=1}^{r}\sigma_{i}({\mathbf{D}}){\mathbf{u}}_{i}{\mathbf{v}}_{i}^{\top}$
be the leading rank-$r$ SVD approximation of ${\mathbf{D}}$ and
${\mathbf{D}}_{\max(-r)}=\sum_{i=r+1}^{m\land
n}\sigma_{i}({\mathbf{D}}){\mathbf{u}}_{i}{\mathbf{v}}_{i}^{\top}$ be its
complement. We also denote ${\rm
SVD}_{r}({\mathbf{D}}):=[{\mathbf{u}}_{1}~{}\cdots{\mathbf{u}}_{r}]$ as the
subspace composed of the leading $r$ left singular vectors of ${\mathbf{D}}$.
The Schatten-$q$ norm of matrix ${\mathbf{D}}$ for $q\geq 1$ is defined as
$\|{\mathbf{D}}\|_{q}:=\left(\sum_{i=1}^{m\land
n}\sigma^{q}_{i}({\mathbf{D}})\right)^{1/q}$. Frobenius norm $\|\cdot\|_{F}$
and spectral norm $\|\cdot\|$ of a matrix are special cases of Schatten-$q$
norm with $q=2$ and $q=\infty$. In addition, ${\mathbf{I}}_{r}$ represents the
$r$-by-$r$ identity matrix. Let
$\mathbb{O}_{p,r}=\\{{\mathbf{U}}:{\mathbf{U}}^{\top}{\mathbf{U}}={\mathbf{I}}_{r}\\}$
be the set of all $p$-by-$r$ matrices with orthonormal columns. For any
${\mathbf{U}}\in\mathbb{O}_{p,r}$,
$P_{{\mathbf{U}}}={\mathbf{U}}{\mathbf{U}}^{\top}$ represents the projection
matrix onto the column span of ${\mathbf{U}}$; we also use
${\mathbf{U}}_{\perp}\in\mathbb{O}_{p,p-r}$ to represent the orthonormal
complement of ${\mathbf{U}}$. We use bracket subscripts to denote sub-
matrices. For example, ${\mathbf{D}}_{[i_{1},i_{2}]}$ is the entry of
${\mathbf{D}}$ on the $i_{1}$-th row and $i_{2}$-th column;
${\mathbf{D}}_{[(r+1):m,:]}$ contains the $(r+1)$-th to the $m$-th rows of
${\mathbf{D}}$. For any matrices ${\mathbf{U}}\in\mathbb{R}^{p_{1}\times
p_{2}}$ and ${\mathbf{V}}\in\mathbb{R}^{m_{1}\times m_{2}}$, let

${\mathbf{U}}\otimes{\mathbf{V}}=\begin{bmatrix}{\mathbf{U}}_{[1,1]}\cdot{\mathbf{V}}&\cdots&{\mathbf{U}}_{[1,p_{2}]}\cdot{\mathbf{V}}\\\
\vdots&&\vdots\\\
{\mathbf{U}}_{[p_{1},1]}\cdot{\mathbf{V}}&\cdots&{\mathbf{U}}_{[p_{1},p_{2}]}\cdot{\mathbf{V}}\\\
\end{bmatrix}\in\mathbb{R}^{(p_{1}m_{1})\times(p_{2}m_{2})}$

be the Kronecker product of ${\mathbf{U}}$ and ${\mathbf{V}}$.

For any order-$d$ tensor
${\mathbfcal{T}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$, let
$\mathcal{M}_{k}(\cdot)$ be the matricization operation that unfolds or
flattens the order-$d$ tensor
${\mathbfcal{T}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$ along mode $k$
into the matrix $\mathcal{M}_{k}({\mathbfcal{T}})\in\mathbb{R}^{p_{k}\times
p_{-k}}$ and here $p_{-k}:=\prod_{j\neq k}p_{j}$. Specifically, the mode-$k$
matricization of ${\mathbfcal{T}}$ is formally defined as

$\mathcal{M}_{k}({\mathbfcal{T}})\in\mathbb{R}^{p_{k}\times
p_{-k}},\quad\left(\mathcal{M}_{k}({\mathbfcal{T}})\right)_{\left[i_{k},j\right]}={\mathbfcal{T}}_{[i_{1},\ldots,i_{d}]},\quad
j=1+\sum_{\begin{subarray}{c}l=1\\\ l\neq
k\end{subarray}}^{d}\left\\{(i_{l}-1)\prod_{\begin{subarray}{c}m=1\\\ m\neq
k\end{subarray}}^{l-1}p_{m}\right\\}$

for any $1\leq i_{l}\leq p_{l},l=1,\ldots,d$. Also see [KB09, Section 2.4] for
more discussion on tensor matricizations. Given two tensors
${\mathbfcal{T}}_{1},{\mathbfcal{T}}_{2}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}}$, define their inner product
$\langle{\mathbfcal{T}}_{1},{\mathbfcal{T}}_{2}\rangle=\sum_{i_{1},\ldots,i_{d}}\langle{\mathbfcal{T}}_{1[i_{1},\ldots,i_{d}]},{\mathbfcal{T}}_{2[i_{1},\ldots,i_{d}]}\rangle$.
The Hilbert-Schmidt norm of ${\mathbfcal{T}}$ is defined as
$\|{\mathbfcal{T}}\|_{{\rm
HS}}=\left(\langle{\mathbfcal{T}},{\mathbfcal{T}}\rangle\right)^{1/2}.$ The
multilinear rank of a tensor ${\mathbfcal{T}}$, ${\rm rank}({\mathbfcal{T}})$,
is defined as a $d$-tuple $(r_{1},\ldots,r_{d})$, where
$r_{k}=\text{rank}(\mathcal{M}_{k}({\mathbfcal{T}}))$ is the mode-$k$ rank.
For any multilinear rank-$(r_{1},\ldots,r_{d})$ tensor ${\mathbfcal{T}}$, it
has Tucker decomposition [Tuc66]:

${\mathbfcal{T}}=\llbracket{\mathbfcal{S}};{\mathbf{U}}_{1},\ldots,{\mathbf{U}}_{d}\rrbracket:={\mathbfcal{S}}\times_{1}{\mathbf{U}}_{1}\times_{2}\cdots\times_{d}{\mathbf{U}}_{d},$
(4)

where ${\mathbfcal{S}}\in\mathbb{R}^{r_{1}\times\cdots\times r_{d}}$ is the
core tensor and ${\mathbf{U}}_{k}\in\mathbb{O}_{p_{i},r_{i}}$ is the mode-$k$
singular subspace. Here, the mode-$k$ product of
${\mathbfcal{T}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$ with a matrix
${\mathbf{U}}\in\mathbb{R}^{p_{k}\times r_{k}}$ is denoted by
${\mathbfcal{T}}\times_{k}{\mathbf{U}}^{\top}$ and is of size
$p_{1}\times\cdots\times p_{k-1}\times r_{k}\times p_{k+1}\times\cdots\times
p_{d}$, such that

$({\mathbfcal{T}}\times_{k}{\mathbf{U}}^{\top})_{[i_{1},\ldots,i_{k-1},j,i_{k+1},\ldots,i_{d}]}=\sum_{i_{k}=1}^{p_{k}}{\mathbfcal{T}}_{[i_{1},i_{2},\ldots,i_{d}]}{\mathbf{U}}_{[i_{k},j]}.$

Given $S=\\{i_{1},\ldots,i_{d}\\}$, it is convenient to denote the product of
${\mathbfcal{T}}$ along the modes indexed by $S$ with the same matrix
${\mathbf{U}}$ and with different matrices $\\{{\mathbf{U}}_{i}\\}$
respectively as

$\begin{split}{\mathbfcal{T}}\times_{S}{\mathbf{U}}:={\mathbfcal{T}}\times_{i_{1}}{\mathbf{U}}\times\cdots\times_{i_{d}}{\mathbf{U}},\quad{\mathbfcal{T}}\times_{i\in
S}{\mathbf{U}}_{i}:={\mathbfcal{T}}\times_{i_{1}}{\mathbf{U}}_{i_{1}}\times\cdots\times_{i_{d}}{\mathbf{U}}_{i_{d}}.\end{split}$

Given
${\mathbfcal{T}}=\llbracket{\mathbfcal{S}};{\mathbf{U}}_{1},\ldots,{\mathbf{U}}_{d}\rrbracket$,
the following relationship of tensor matricization is used frequently in the
proof:

$\mathcal{M}_{k}\left({\mathbfcal{S}}\times_{1}{\mathbf{U}}_{1}\times_{2}\cdots\times_{d}{\mathbf{U}}_{d}\right)={\mathbf{U}}_{k}\mathcal{M}_{k}({\mathbfcal{S}})({\mathbf{U}}_{d}^{\top}\otimes\cdots\otimes{\mathbf{U}}_{k+1}^{\top}\otimes{\mathbf{U}}_{k-1}^{\top}\otimes\cdots\otimes{\mathbf{U}}_{1}^{\top}).$
(5)

We refer the readers to [Kol01, KB09] for a more comprehensive survey on
tensor algebra.

Finally, we use $\sin\Theta$ distance to measure the difference between two
$p$-by-$r$ column orthogonal matrices $\widehat{{\mathbf{U}}}$ and
${\mathbf{U}}$. Suppose the singular values of
$\widehat{{\mathbf{U}}}^{\top}{\mathbf{U}}$ are
$\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{r}\geq 0$. Then
$\Theta(\widehat{{\mathbf{U}}},{\mathbf{U}})$ is defined as

$\Theta(\widehat{{\mathbf{U}}},{\mathbf{U}})={\rm
diag}\left(\cos^{-1}(\sigma_{1}),\cos^{-1}(\sigma_{2}),\ldots,\cos^{-1}(\sigma_{r})\right).$

Common properties of $\sin\Theta$ distance can be found in [CZ18, Lemma 1] and
[LZ20a, Lemma 7].

###  2.1 Blockwise Errors of ${\mathbfcal{Z}}$

In this subsection, we introduce key quantities appearing in the perturbation
bounds that characterize the blockwise errors of ${\mathbfcal{Z}}$. For
simplicity, we consider order-3 tensors and
$\Omega_{1}=\\{1\\},\Omega_{2}=\\{2\\},\Omega_{3}=\\{3\\}$ perturbation
setting for illustration.

Define the blockwise errors of ${\mathbfcal{Z}}$ that characterize the tensor
perturbation:

$\begin{split}&\tau_{1}=\max_{k\in[3]}\tau_{1k},\quad\tau_{1k}=\left\|\left(\mathcal{M}_{k}({\mathbfcal{Z}}\times_{k+1}{\mathbf{U}}_{k+1}^{\top}\times_{k+2}{\mathbf{U}}_{k+2}^{\top})\right)_{\max(r_{k})}\right\|_{q},\quad
k=1,2,3;\\\
&\tau_{2}=\max_{k\in[3]}\Big{\\{}\max_{\begin{subarray}{c}{\mathbf{V}}\in\mathbb{R}^{(p_{k+1}-r_{k+1})\times
r_{k+1}}\\\ \|{\mathbf{V}}\|_{q}\leq
1\end{subarray}}\left\|\left(\mathcal{M}_{k}({\mathbfcal{Z}}\times_{k+1}({\mathbf{U}}_{k+1\perp}{\mathbf{V}})^{\top}\times_{k+2}{\mathbf{U}}_{k+2}^{\top})\right)_{\max(r_{k})}\right\|_{q},\\\
&\quad\quad\quad\quad\max_{\begin{subarray}{c}{\mathbf{V}}\in\mathbb{R}^{(p_{k+2}-r_{k+2})\times
r_{k+2}}\\\ \|{\mathbf{V}}\|_{q}\leq
1\end{subarray}}\left\|\left(\mathcal{M}_{k}({\mathbfcal{Z}}\times_{k+1}{\mathbf{U}}_{k+1}^{\top}\times_{k+2}({\mathbf{U}}_{k+2\perp}{\mathbf{V}})^{\top})\right)_{\max(r_{k})}\right\|_{q}\Big{\\}};\\\
&\tau_{3}=\max_{k\in[3]}\max_{\begin{subarray}{c}{\mathbf{V}}\in\mathbb{R}^{(p_{k+1}-r_{k+1})\times
r_{k+1}}\\\ {\mathbf{V}}^{\prime}\in\mathbb{R}^{(p_{k+2}-r_{k+2})\times
r_{k+2}}\\\ \|{\mathbf{V}}\|_{q}\leq 1,\|{\mathbf{V}}^{\prime}\|_{q}\leq
1\end{subarray}}\left\|\left(\mathcal{M}_{k}({\mathbfcal{Z}}\times_{k+1}({\mathbf{U}}_{k+1\perp}{\mathbf{V}})^{\top}\times_{k+2}({\mathbf{U}}_{k+2\perp}{\mathbf{V}}^{\prime})^{\top})\right)_{\max(r_{k})}\right\|_{q}.\end{split}$
(6)

Here all mode indices $(\cdot)_{k}$ of an order-3 tensor are in the sense of
modulo-3, e.g., $r_{1}=r_{4}$, $p_{2}=p_{5}$.

$\tau_{1},\tau_{2},\tau_{3}$ represent the maximum of blockwise errors of
${\mathbfcal{Z}}$ in the projection spaces expanded by
${\mathbf{U}}_{1},{\mathbf{U}}_{2},{\mathbf{U}}_{3}$ and their complements.
For example, in Figure 1 we illustrate the blockwise errors characterized by
$\tau_{11},\tau_{12},\tau_{13}$. $\tau_{2},\tau_{3}$ characterize blockwise
errors of ${\mathbfcal{Z}}$ in a similar fashion but with more complicated
projections.

(a) $\tau_{11}$ (b) $\tau_{12}$ (c) $\tau_{13}$

Figure 1: Illustration of $\tau_{11},\tau_{12},\tau_{13}$. Here, we assume
${\mathbf{U}}_{k}^{\top}=[{\mathbf{I}}_{r_{k}}~{}\boldsymbol{0}_{r_{k}\times(p_{k}-r_{k})}]$,
$k=1,2,3$, for a better visualization. The red, green, blue blocks represent
the corresponding blockwise errors $\tau_{11},\tau_{12},\tau_{13}$ in
${\mathbfcal{Z}}$.

These blockwise errors of ${\mathbfcal{Z}}$ are in fact a generalization of
error terms in matrix perturbation theory. In the matrix setting
$\widetilde{{\mathbf{T}}}={\mathbf{T}}+{\mathbf{Z}}$, let
$\widehat{{\mathbf{U}}},\widehat{{\mathbf{V}}}$ and
${\mathbf{U}},{\mathbf{V}}$ be leading left and right singular vectors of
$\widetilde{{\mathbf{T}}}$ and ${\mathbf{T}}$, respectively. Then by Wedin’s
perturbation theory, the upper bounds of
$\left\|\sin\Theta(\widehat{{\mathbf{U}}},{\mathbf{U}})\right\|$ and
$\left\|\sin\Theta(\widehat{{\mathbf{V}}},{\mathbf{V}})\right\|$ involve
$\|{\mathbf{Z}}\widehat{{\mathbf{V}}}\|$ and
$\|\widehat{{\mathbf{U}}}^{\top}{\mathbf{Z}}\|$, which are also blockwise
errors of ${\mathbf{Z}}$.

Next, we introduce a simple quantity $\xi$ that characterizes the error bound
for tensor reconstruction. In this order-3 asymmetric setting, $\xi$ is
defined as

$\xi:=\underset{\|{\mathbfcal{Y}}\|_{{\rm HS}}\leq 1,{\rm
rank}({\mathbfcal{Y}})\leq(r_{1},r_{2},r_{3})}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\rangle.$
(7)

In the following Lemma 1, we give several equivalent characterizations of
$\xi$.

######  Lemma 1 (Equivalent Characterizations of $\xi$)

$\begin{split}\xi&:=\underset{\|{\mathbfcal{Y}}\|_{{\rm HS}}\leq 1,{\rm
rank}({\mathbfcal{Y}})\leq(r_{1},r_{2},r_{3})}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\rangle\\\
&=\sup_{{\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}},1\leq i\leq
3}\|{\mathbfcal{Z}}\times_{1}{\mathbf{U}}_{1}^{\top}\times{\mathbf{U}}_{2}^{\top}\times_{d}{\mathbf{U}}_{3}^{\top}\|_{{\rm
HS}}\\\ &=\inf_{{\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}},1\leq i\leq
3}\|{\mathbfcal{Z}}-{\mathbfcal{Z}}\times_{1}P_{{\mathbf{U}}_{1}}\times
P_{{\mathbf{U}}_{2}}\times_{d}P_{{\mathbf{U}}_{3}}\|_{{\rm HS}}.\end{split}$

By Lemma 1, $\xi$ measures the residual between ${\mathbfcal{Z}}$ and its best
multilinear rank-$(r_{1},r_{2},r_{3})$ approximation in Hilbert-Schmidt norm.
Another interpretation of $\xi$ is from Gaussian width [Gor88], which will be
discussed in Section 6.

Although the exact values of $\tau_{j}$ and $\xi$ may be hard to compute in
general, it is often practical to provide probabilistic bounds when we impose
distributional assumptions on ${\mathbfcal{Z}}$. For example, if
${\mathbfcal{Z}}$ is a random tensor with i.i.d. standard normal entries and
consider $r_{1}=r_{2}=r_{3}=r,p_{1}=p_{2}=p_{3}=p$, $q=\infty$, then by random
matrix theory [Ver10], we can show that with high probability
$\tau_{1i}\asymp(\sqrt{p}+r)$,
$\tau_{2}\asymp\sqrt{pr},\tau_{3}\asymp\sqrt{pr}$ and
$\xi\asymp(r^{\frac{3}{2}}+\sqrt{pr})$. More details about the upper bounds
for $\tau_{j}$ and $\xi$ can be found in tensor denoising and tensor co-
clustering applications in Section 6.

##  3 Illustration of Perturbation Bounds for $d=3$ Asymmetric Case

In this section, we present our main results in the $d=3$ asymmetric case to
better illustrate the main ideas in this paper. The specialized HOOI algorithm
for the $d=3$ asymmetric case is given in Algorithm 2 and its guarantee is
provided in Theorem 1.

Algorithm 2 High-Order Orthogonal Iteration for Tensor Decomposition ($d=3$)

Input: $\widetilde{{\mathbfcal{T}}}\in\mathbb{R}^{p_{1}\times p_{2}\times
p_{3}}$, initialization $\\{\widetilde{{\mathbf{U}}}_{i}^{(0)}\\}_{i=1}^{3}$
with $\widetilde{{\mathbf{U}}}_{i}^{(0)}\in\mathbb{O}_{p_{i},r_{i}}$, maximum
number of iterations $t_{\max}$.
Output:
$\\{\widehat{{\mathbf{U}}}_{i}\\}_{i=1}^{3},\widehat{{\mathbfcal{T}}}$.

1:For $t=1,\ldots,t_{\max}$, do
$\begin{split}&\widetilde{{\mathbf{U}}}_{1}^{(t+1)}={\rm
SVD}_{r_{1}}\left(\mathcal{M}_{1}\left(\widetilde{{\mathbfcal{T}}}\times_{2}(\widetilde{{\mathbf{U}}}_{2}^{(t)})^{\top}\times_{3}(\widetilde{{\mathbf{U}}}_{3}^{(t)})^{\top}\right)\right),\\\
&\widetilde{{\mathbf{U}}}_{2}^{(t+1)}={\rm
SVD}_{r_{2}}\left(\mathcal{M}_{2}\left(\widetilde{{\mathbfcal{T}}}\times_{1}(\widetilde{{\mathbf{U}}}_{1}^{(t+1)})^{\top}\times_{3}(\widetilde{{\mathbf{U}}}_{3}^{(t)})^{\top}\right)\right),\\\
&\widetilde{{\mathbf{U}}}_{3}^{(t+1)}={\rm
SVD}_{r_{3}}\left(\mathcal{M}_{3}\left(\widetilde{{\mathbfcal{T}}}\times_{1}(\widetilde{{\mathbf{U}}}_{1}^{(t+1)})^{\top}\times_{2}(\widetilde{{\mathbf{U}}}_{2}^{(t+1)})^{\top}\right)\right).\end{split}$

2:Let $\widehat{{\mathbf{U}}}_{i}=\widetilde{{\mathbf{U}}}_{i}^{(t_{\max})}$
for $i=1,2,3$ and compute
$\widehat{{\mathbfcal{T}}}=\widetilde{{\mathbfcal{T}}}\times_{1}P_{\widehat{{\mathbf{U}}}_{1}}\times_{2}P_{\widehat{{\mathbf{U}}}_{2}}\times_{3}P_{\widehat{{\mathbf{U}}}_{3}}.$

######  Theorem 1 (Tensor Perturbation Bounds for HOOI (d = 3))

Consider the perturbation model (1) with
$\widetilde{{\mathbfcal{T}}},{\mathbfcal{T}},{\mathbfcal{Z}}\in\mathbb{R}^{p_{1}\times
p_{2}\times p_{3}}$. Define the blockwise errors as in (6) (7) and denote the
initialization error of $\\{\widetilde{{\mathbf{U}}}_{k}^{(0)}\\}_{k=1}^{3}$
as
$e_{0}:=\max_{k=1,2,3}\|\widetilde{{\mathbf{U}}}^{(0)\top}_{k\perp}{\mathbf{U}}_{k}\|_{q}$.
Assume the initialization error and the signal strength satisfy

$e_{0}\leq\sqrt{2}/2\,\text{ and }\,\lambda\geq(16+12\sqrt{2})\xi.$ (8)

Let
$\widetilde{{\mathbfcal{T}}}^{(t)}:=\widetilde{{\mathbfcal{T}}}\times_{1}P_{\widetilde{{\mathbf{U}}}^{(t)}_{1}}\times_{2}P_{\widetilde{{\mathbf{U}}}^{(t)}_{2}}\times_{3}P_{\widetilde{{\mathbf{U}}}^{(t)}_{3}}$
be the estimator of ${\mathbfcal{T}}$ after $t$ steps in Algorithm 1. Then
with inputs $\widetilde{{\mathbfcal{T}}}$,
$\\{\widetilde{{\mathbf{U}}}_{i}^{(0)}\\}_{i=1}^{3}$, the mode-$k$ singular
subspace updates in Algorithm 2 after $t$ iterations satisfy

$\max_{k=1,2,3}\left\|\sin\Theta\left(\widetilde{{\mathbf{U}}}^{(t)}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}\leq\frac{8\tau_{1}}{\lambda}+\frac{e_{0}}{2^{t}},$
(9)

and the $t$-step tensor estimation satisfies

$\left\|\widetilde{{\mathbfcal{T}}}^{(t)}-{\mathbfcal{T}}\right\|_{{\rm
HS}}\leq\left(1+\frac{6}{1-\left(8\frac{\tau_{1}}{\lambda}+\frac{e_{0}}{2^{t-1}}\right)^{2}}\right)\xi.$
(10)

Moreover, when $t_{\max}\geq C\log(e_{0}\lambda/\tau_{1})\vee 1$ for some
$C>0$, the outputs of the estimated mode-$k$ singular subspaces of Algorithm 2
satisfy

$\max_{k=1,2,3}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}\leq\frac{9\tau_{1}}{\lambda},$
$\begin{split}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}=\left\|\widehat{{\mathbf{U}}}_{k\perp}^{\top}{\mathbf{U}}_{k}\right\|_{q}\leq
4\left(\frac{\tau_{1k}}{\lambda}+\frac{18\tau_{1}\tau_{2}}{\lambda^{2}}+\frac{81\tau_{1}^{2}\tau_{3}}{\lambda^{3}}\right),\quad
k=1,2,3,\end{split}$ (11)

and the output of tensor reconstruction $\widehat{{\mathbfcal{T}}}$ satisfies

$\begin{split}\left\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\right\|_{{\rm
HS}}&\leq\left\|\llbracket{\mathbfcal{Z}};\widehat{{\mathbf{U}}}_{1}^{\top},\widehat{{\mathbf{U}}}_{2}^{\top},\widehat{{\mathbf{U}}}_{3}^{\top}\rrbracket\right\|_{{\rm
HS}}+\sum_{k=1}^{3}\left\|\widehat{{\mathbf{U}}}_{k\perp}^{\top}\mathcal{M}_{k}({\mathbfcal{T}})\right\|_{F}\leq
13\xi.\end{split}$ (12)

######  Remark 1 (Initialization Error $e_{0}$ and Signal Strength $\lambda$)

In Theorem 1, we assume the initializations
$\\{\widetilde{{\mathbf{U}}}_{i}^{(0)}\\}_{i=1}^{3}$ are adequate, i.e., the
maximum error $e_{0}$ is bounded by a constant. Here $\frac{\sqrt{2}}{2}$ is
chosen for convenience and the same proof works for any fixed constant less
than $1$. The initialization can be computed by using different methods
depending on the application. For example, in tensor PCA/SVD [RM14, ZX18], we
can let $\widetilde{{\mathbf{U}}}_{i}^{(0)}={\rm
SVD}_{r_{i}}(\mathcal{M}_{i}(\widetilde{{\mathbfcal{T}}}))$; in tensor
completion [XYZ17], we can choose $\widetilde{{\mathbf{U}}}_{i}^{(0)}$ to be
the leading $r_{i}$ singular vectors of
$\mathcal{M}_{i}(\widetilde{{\mathbfcal{T}}})\mathcal{M}_{i}(\widetilde{{\mathbfcal{T}}})^{\top}$
with diagonal deletion.

The signal strength requirement on $\lambda$ is in the same vein as the
singular gap condition appearing in the classical matrix perturbation theory.

######  Remark 2 (Mode-$k$ Singular Subspace Linear Convergence Property)

From (9), we know after $t$ steps, the upper bound of
$e_{t}:=\max_{k=1,2,3}\left\|\sin\Theta\left(\widetilde{{\mathbf{U}}}^{(t)}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}$
can be decomposed into two parts. The first part is a fixed quantity over
iterations which represents the intrinsic estimation error due to the
perturbation of the observed data. The second part decreases linearly to $0$
with respect to the iteration number and could be viewed as the avoidable
error by leveraging the low multilinear rank structure of ${\mathbfcal{T}}$.

######  Remark 3 (Unilateral Perturbation Bounds for Mode-$k$ Singular
Subspace)

Our tensor perturbation bounds on singular subspace share the same spirit as
the unilateral perturbation bounds on singular subspace of matrix SVD in
[CZ18]. Consider the matrix perturbation setting mentioned in Section 2.1 with
the additional assumption that ${\mathbf{T}}$ is rank-$r$ and has SVD
${\mathbf{U}}\boldsymbol{\Sigma}{\mathbf{V}}^{\top}$. [CZ18] showed that the
upper bound of $\|\sin\Theta(\widehat{{\mathbf{U}}},{\mathbf{U}})\|$ can be
written as $\frac{a_{1}}{\lambda}+\frac{a_{2}}{\lambda^{2}}$, which can be
interpreted as the sum of first and second order perturbations. In Theorem 1,
the upper bound of
$\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}$
can be also written as $\sum_{i=1}^{3}\frac{a_{i}}{\lambda^{i}}$ which can be
interpreted as summation of the first, second, and third order perturbations.
This phenomenon also generalizes to order-$d$ case in Theorem 3.

Due to the unilateral property, when the tensor dimension of each mode is at
different order, the estimation error rate of singular subspace in each mode
can vary significantly. For example in the tensor denoising setting,
$\widetilde{{\mathbfcal{T}}}={\mathbfcal{T}}+{\mathbfcal{Z}}$ where
${\mathbfcal{T}}\in\mathbb{R}^{p_{1}\times p_{2}\times p_{3}}$ is a
multilinear rank-$(r_{1},r_{2},r_{3})$ tensor and ${\mathbfcal{Z}}$ is a
random tensor with i.i.d. standard normal entries. Let
$r_{\max}=\max_{i}r_{i}$ and suppose $p_{1}\ll p_{2}\ll p_{3}$, $r_{\max}\ll
p_{1}^{1/2}$. Consider $q=\infty$, then by random matrix theory [Ver10], we
can show $\tau_{1i}\leq C\sqrt{p_{i}}$, $\tau_{2}\leq C(\sqrt{p_{3}r_{\max}})$
and $\tau_{3}\leq C(\sqrt{p_{3}r_{\max}})$ with high probability. Thus when
$\lambda\geq Cp_{3}\sqrt{\frac{r_{\max}}{p_{1}}}$, Theorem 1 immediately
implies, with high probability

$\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|\leq
C\frac{\sqrt{p_{k}}}{\lambda},\quad k=1,2,3$

for some $C>0$. So we can see the perturbation of $\widehat{{\mathbf{U}}}_{k}$
depends on $p_{k}$. Also as $\lambda$ decreases, for different $k$, different
order perturbations in (11) could dominate in the perturbation bound of
$\widehat{{\mathbf{U}}}_{k}$. For example, when $\lambda\asymp
Cp_{3}\sqrt{\frac{r_{\max}}{p_{2}}}$, Theorem 1 yields

$\begin{split}&\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{1},{\mathbf{U}}_{1}\right)\right\|\leq
C\frac{\sqrt{p_{3}^{2}r_{\max}}}{\lambda^{2}},\\\
&\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{2},{\mathbf{U}}_{2}\right)\right\|\leq
C\frac{\sqrt{p_{2}}}{\lambda}+C^{\prime}\frac{\sqrt{p_{3}^{2}r_{\max}}}{\lambda^{2}},\\\
&\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{3},{\mathbf{U}}_{3}\right)\right\|\leq
C\frac{\sqrt{p_{3}}}{\lambda}\end{split}$

for constants $C,C^{\prime}>0$. More details about the application of HOOI
perturbation bounds in order-$d$ tensor denoising and numerical studies for
this unilateral property of singular subspace perturbation can be found in
Sections 6.1 and 7.1, respectively.

######  Remark 4 (Comparison of Perturbation Bounds of truncated HOSVD and
HOOI)

It is worth mentioning that the power iteration in Algorithm 1 plays an
important role for refining tensor reconstruction. Without power iteration,
the estimator $\widehat{{\mathbfcal{T}}}^{{\rm
T-HOSVD}}=\widetilde{{\mathbfcal{T}}}\times_{1}P_{\widetilde{{\mathbf{U}}}_{1}^{(0)}}\times_{2}P_{\widetilde{{\mathbf{U}}}_{2}^{(0)}}\times_{3}P_{\widetilde{{\mathbf{U}}}_{3}^{(0)}}$
with $\widetilde{{\mathbf{U}}}_{i}^{(0)}={\rm
SVD}_{r_{i}}(\mathcal{M}_{i}(\widetilde{{\mathbfcal{T}}}))$ is called
truncated HOSVD (T-HOSVD) in the literature [DLDMV00a]. It is not hard to show
$\|\widehat{{\mathbfcal{T}}}^{{\rm T-HOSVD}}-{\mathbfcal{T}}\|_{{\rm HS}}\leq
C\|{\mathbfcal{Z}}\|_{{\rm HS}}$ for some $C>0$. Since
$\|{\mathbfcal{Z}}\|_{{\rm HS}}=\underset{\|{\mathbfcal{Y}}\|_{{\rm HS}}\leq
1}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\rangle$ and may be much larger
than $\xi$, we can see the power iteration can greatly improve the accuracy
for tensor reconstruction, and this echos the findings in literature in tensor
denoising setting [ZX18].

The following lemma provides an alternative way to bound
$\left\|\llbracket{\mathbfcal{Z}};\widehat{{\mathbf{U}}}_{1}^{\top},\widehat{{\mathbf{U}}}_{2}^{\top},\widehat{{\mathbf{U}}}_{3}^{\top}\rrbracket\right\|_{{\rm
HS}}$ appearing in the reconstruction error bound (12).

###### Lemma 2

Suppose ${\mathbfcal{Z}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$ is a
general order-$d$ tensor and
${\mathbf{U}}_{k},\widehat{{\mathbf{U}}}_{k}\in\mathbb{O}_{p_{k},r_{k}}$ are
general matrices with orthonormal columns. For any subset
$\Omega\subseteq\\{1,\ldots,d\\}$, we further define projections of
${\mathbfcal{Z}}$ on $\Omega$ as follows,

$\theta_{\Omega}=\left\|{\mathbfcal{Z}}\times_{k\in\Omega}{\mathbf{U}}_{k}^{\top}\times_{k\in\Omega^{c}}{\mathbf{U}}_{k\perp}^{\top}\right\|_{{\rm
HS}}.$

Then,

$\left\|\llbracket{\mathbfcal{Z}};\widehat{{\mathbf{U}}}_{1}^{\top},\widehat{{\mathbf{U}}}_{2}^{\top},\ldots,\widehat{{\mathbf{U}}}_{d}^{\top}\rrbracket\right\|_{{\rm
HS}}\leq\sum_{\Omega\subseteq\\{1,\ldots,d\\}}\theta_{\Omega}\prod_{k\in\Omega^{c}}\left\|\sin\Theta(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k})\right\|.$

We end this section by introducing a deterministic rate matching lower bound
for tensor reconstruction. Since the statement of the lower bound is relative
simple, we state it in general order-$d$ setting. In particular, we consider
the following class of $({\mathbfcal{T}},{\mathbfcal{Z}})$ pairs of
$p_{1}\times\cdots\times p_{d}$ tensors and perturbations,

$\mathcal{F}_{r}(\xi)=\left\\{\left({\mathbfcal{T}},{\mathbfcal{Z}}\right):{\rm
rank}({\mathbfcal{T}})=(r_{1},\cdots,r_{d})\leq{\mathbf{r}},\underset{\begin{subarray}{c}\|{\mathbfcal{Y}}\|_{{\rm
HS}}\leq 1,\\\ {\rm
rank}({\mathbfcal{Y}})\leq(r_{1},\ldots,r_{d})\end{subarray}}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\rangle\leq\xi\right\\},$

here ${\mathbf{r}}=(r,\ldots,r)$ and the comparison is entrywise.

######  Theorem 2 (Tensor Reconstruction Lower Bound under Perturbation)

Consider perturbation model (1), we have the following deterministic lower
bound for reconstructing ${\mathbfcal{T}}$,

$\inf_{\widehat{{\mathbfcal{T}}}}\sup_{({\mathbfcal{T}},{\mathbfcal{Z}})\in\mathcal{F}_{r}(\xi)}\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm
HS}}\geq\frac{\sqrt{2}}{2}\xi.$

######  Remark 5 (Optimality of HOOI and one-step HOOI for Tensor
Reconstruction)

When tensor order $d$ is fixed, combining Theorem 1 and 2, we have shown that
HOOI with good initialization is optimal for tensor reconstruction in the
class $({\mathbfcal{T}},{\mathbfcal{Z}})\in\mathcal{F}_{r}(\xi)$. At the same
time, from (10), we see the error rate of tensor reconstruction is optimal
even after one iteration of HOOI i.e., $t_{max}=1$ and more iterations can
improve the coefficient in front of $\xi$. This suggests that in some
applications where running HOOI until convergence is prohibitive, we can just
run it for one iteration to get a fairly good reconstruction. See more in
Section 7.2 about a numerical comparison of HOOI and one-step HOOI.

##  4 A Blockwise Perturbation Bound of Higher-order Orthogonal Iteration for
Tensor Decomposition

In this section, we present the main results of perturbation bounds of HOOI
given in Algorithm 1. In contrast with Theorem 1, Theorem 3 in this section
covers the general order-$d$ perturbation setting with
$\widetilde{{\mathbfcal{T}}}$ having symmetric index groups
$\\{\Omega_{i}\\}_{i=1}^{m}$. Before stating the theorem, we first define the
blockwise errors of ${\mathbfcal{Z}}$ in this general setting. Let
$\mathcal{S}_{i}^{(-\bar{k})}:=\left\\{S\subseteq[d]\setminus\\{\bar{k}\\}:|S|=i\right\\}$
be the collection of all possible index sets with $i$ elements from
$[d]\setminus\\{\bar{k}\\}$ and $\mathcal{S}_{0}^{(-\bar{k})}:=\emptyset$. For
$S\in\mathcal{S}_{i}^{(-\bar{k})}$, we let
$S^{c}=([d]\setminus\\{\bar{k}\\})\setminus S$. Now we define the blockwise
errors of ${\mathbfcal{Z}}$ as

$\begin{split}&\tau_{1}=\max_{k=1,\ldots,m}\tau_{1k},\quad\tau_{1k}=\left\|\left(\mathcal{M}_{\bar{k}}({\mathbfcal{Z}}\times_{i\neq\bar{k}}{\mathbf{U}}_{i^{\prime}}^{\top})\right)_{\max(r_{k})}\right\|_{q},\quad
k=1,\ldots,m;\\\
&\tau_{j}=\max_{k\in[m]}\Big{\\{}\max_{S\in\mathcal{S}^{(-\bar{k})}_{j-1}}\sup_{\begin{subarray}{c}{\mathbf{V}}_{i^{\prime}}\in\mathbb{R}^{(p_{i^{\prime}}-r_{i^{\prime}})\times
r_{i^{\prime}}},\\\ \|{\mathbf{V}}_{i^{\prime}}\|_{q}\leq 1,i\in
S\end{subarray}}\left\|\left(\mathcal{M}_{\bar{k}}({\mathbfcal{Z}}\times_{i\in
S}\left({\mathbf{U}}_{i^{\prime}\perp}{\mathbf{V}}_{i^{\prime}}\right)^{\top}\times_{i\in
S^{c}}{\mathbf{U}}_{i^{\prime}}^{\top}\right)_{\max(r_{k})}\right\|_{q}\Big{\\}},\\\
&\quad\text{ for }j=2,\ldots,m.\end{split}$ (13)

Finally, $\xi$ in this setting is defined as

$\xi:=\underset{\|{\mathbfcal{Y}}\|_{{\rm HS}}\leq 1,{\rm
rank}({\mathbfcal{Y}})\leq(r_{1^{\prime}},\ldots,r_{d^{\prime}})}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\rangle.$
(14)

######  Theorem 3 (General Perturbation Bounds for Tensor Power Iteration)

Consider the perturbation model (1) with
$\widetilde{{\mathbfcal{T}}},{\mathbfcal{T}},{\mathbfcal{Z}}\in\mathbb{R}^{p_{1^{\prime}}\times\cdots\times
p_{d^{\prime}}}$, symmetric index groups $(\Omega_{1},\ldots,\Omega_{m})$ and
blockwise errors in (13) (14). Denote the initialization error of
$\\{\widetilde{{\mathbf{U}}}_{k}^{(0)}\\}_{k=1}^{m}$ as
$e_{0}:=\max_{k=1,\ldots,m}\|\widetilde{{\mathbf{U}}}^{(0)\top}_{k\perp}{\mathbf{U}}_{k}\|_{q}$.
Assume the initialization error and the signal strength satisfy

$e_{0}\leq\frac{\sqrt{2}}{2}\,\text{ and }\,\lambda\geq
2^{\frac{d+4}{2}}\left(1+\frac{\sqrt{2}}{2}\right)^{d}\xi.$ (15)

Let
$\widetilde{{\mathbfcal{T}}}^{(t)}:=\widetilde{{\mathbfcal{T}}}\times_{\Omega_{1}}P_{\widetilde{{\mathbf{U}}}^{(t)}_{1}}\times\cdots\times_{\Omega_{m}}P_{\widetilde{{\mathbf{U}}}^{(t)}_{m}}$
be the estimator of ${\mathbfcal{T}}$ after $t$ steps in Algorithm 1. Then
with inputs $\widetilde{{\mathbfcal{T}}}$,
$\\{\widetilde{{\mathbf{U}}}_{k}^{(0)}\\}_{k=1}^{m}$,
$\\{\Omega_{i}\\}_{i=1}^{m}$, the mode-$k$ singular subspace updates in
Algorithm 1 after $t$ iterations satisfy

$\max_{k\in[m]}\left\|\sin\Theta\left(\widetilde{{\mathbf{U}}}_{k}^{(t)},{\mathbf{U}}_{k}\right)\right\|_{q}\leq
2^{\frac{d+3}{2}}\frac{\tau_{1}}{\lambda}+\frac{e_{0}}{2^{t}},$ (16)

and the $t$-step tensor estimation satisfies

$\left\|\widetilde{{\mathbfcal{T}}}^{(t)}-{\mathbfcal{T}}\right\|_{{\rm
HS}}\leq\left(1+2d\left(1-\left(2^{\frac{d+3}{2}}\frac{\tau_{1}}{\lambda}+\frac{e_{0}}{2^{t-1}}\right)^{2}\right)^{-\frac{d-1}{2}}\right)\xi.$
(17)

Moreover, when $t_{\max}\geq C\log(e_{0}\lambda/\tau_{1})\vee 1$ for some
$C>0$, the outputs of estimated mode-$k$ singular subspace of Algorithm 1
satisfy

$\max_{k\in[m]}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}\leq\left(2^{\frac{d+3}{2}}+1\right)\frac{\tau_{1}}{\lambda},$
$\begin{split}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}&\leq\frac{2}{\left(1-c^{*}(\tau_{1},\lambda,d)\right)^{\frac{d-1}{2}}}\left(\frac{\tau_{1k}}{\lambda}+\sum_{j=1}^{d-1}\frac{{d-1\choose
j}\left(2^{\frac{d+3}{2}}+1\right)^{j}\tau_{1}^{j}\tau_{j+1}}{\lambda^{j+1}}\right).\end{split}$
(18)

for $k=1,\ldots,m$ where
$c^{*}(\tau_{1},\lambda,d):=\frac{\left(2^{\frac{d+3}{2}}+2\right)^{2}\tau_{1}^{2}}{\lambda^{2}}\leq\frac{1}{2}$,
and the output of tensor reconstruction $\widehat{{\mathbfcal{T}}}$ satisfies

$\begin{split}\left\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\right\|_{{\rm
HS}}\leq&\left\|{\mathbfcal{Z}}\times_{\Omega_{1}}\widehat{{\mathbf{U}}}_{1}^{\top}\times\cdots\times_{\Omega_{m}}\widehat{{\mathbf{U}}}_{m}^{\top}\right\|_{{\rm
HS}}+\sum_{k=1}^{m}|\Omega_{k}|\left\|\widehat{{\mathbf{U}}}_{k\perp}^{\top}\mathcal{M}_{\bar{k}}({\mathbfcal{T}})\right\|_{F}\\\
\leq&\left(1+2d\left(1-c^{*}\left(\tau_{1},\lambda,d\right)\right)^{-\frac{d-1}{2}}\right)\xi.\end{split}$
(19)

######  Remark 6 (Size of $c^{*}(\tau_{1},\lambda,d)$)

It is easy to check $c^{*}(\tau_{1},\lambda,d)\leq\frac{1}{2}$ based on
$\tau_{1}\leq\xi$ and the requirement of the signal strength $\lambda$. So we
have $\left(1-c^{*}(\tau_{1},\lambda,d)\right)^{-\frac{d-1}{2}}\leq
2^{\frac{d-1}{2}}$ in the upper bounds of
$\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}$
and $\left\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\right\|_{{\rm HS}}$.
However, in many practical applications, such as tensor denoising to be
introduced in Section 6.1, tensor order $d$ is fixed and
$\lambda\gg(2^{\frac{d+3}{2}}+2)\tau_{1}$. In this case
$c^{*}(\tau_{1},\lambda,d)$ could be much smaller than $\frac{1}{2}$ and the
scale of $\left(1-c^{*}(\tau_{1},\lambda,d)\right)^{-\frac{d-1}{2}}$ can be
very close to $1$.

######  Remark 7 (Comparison with [AGH+14])

Compared with the perturbation bounds of power iteration for supersymmetric
CP-low-rank decomposition [AGH+14, Theorem 5.1], our Theorem 3 covers more
general symmetric and partial symmetric multilinear low-rank decomposition
settings. Also in Theorem 5.1 of [AGH+14], the tensor reconstruction error
bound of power iteration is given in terms of tensor spectral norm, which does
not improve upon the guarantee by the trivial estimator
$\widetilde{{\mathbfcal{T}}}$. On the other hand, the tensor reconstruction
error of $\widehat{{\mathbfcal{T}}}$ in Theorem 3 is given in Hilbert-Schmidt
norm and can be significantly better than the guarantee for
$\widetilde{{\mathbfcal{T}}}$ as $\|{\mathbfcal{Z}}\|_{{\rm HS}}\gg\xi$ in
most of the applications.

######  Remark 8 (A Proof Sketch of Theorem 3)

We provide a sketch on how to prove (16) and (17). The rest of the results
(18) and (19) follow easily from (16) (17) by plugging in
$t_{\max}\geq\log(e_{0}\lambda/\tau_{1})\vee 1$. The idea is to develop the
recursive error bounds of $\widetilde{{\mathbf{U}}}_{k}^{(t+1)}$ i.e., the
estimate of ${\mathbf{U}}_{k}$ at iteration $t+1$, based on the error bound of
$\widetilde{{\mathbf{U}}}_{k}^{(t)}$, i.e., the estimate at iteration $t$. The
argument can be divided into three steps. It is worth mentioning that all
three steps involves complex tensor algebra and this makes the proof even more
difficult.

First, we denote

$e_{t}=\max_{k}e_{t,k},\quad
e_{t,k}=\left\|(\widetilde{{\mathbf{U}}}_{k\perp}^{(t)})^{\top}{\mathbf{U}}_{k}\right\|_{q},\quad
k=1,\ldots,m;t=0,1,\ldots.$

Step 1. In HOOI procedure, the update for the mode-$k$ singular subspace
satisfies

$\begin{split}\widetilde{{\mathbf{U}}}_{k}^{(t+1)}=&{\rm
SVD}_{r_{k}}\bigg{(}\mathcal{M}_{\bar{k}}\left({\mathbfcal{T}}\times_{i\in\underline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t+1)\top}_{i^{\prime}}\times_{\widecheck{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{k}\times_{i\in\overline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{i^{\prime}}\right)\\\
&\quad\quad\quad+\mathcal{M}_{\bar{k}}\left({\mathbfcal{Z}}\times_{i\in\underline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t+1)\top}_{i^{\prime}}\times_{\widecheck{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{k}\times_{i\in\overline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{i^{\prime}}\right)\bigg{)},\end{split}$

here
$\underline{\Omega}_{i}:=\bigcup_{j=1}^{i-1}\Omega_{j},\overline{\Omega}_{i}:=\bigcup_{j=i+1}^{d}\Omega_{j}$.
To give an upper bound for $e_{t+1,k}$, we aim to give an upper bound for

$\left\|\left(\mathcal{M}_{\bar{k}}\left({\mathbfcal{Z}}\times_{i\in\underline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t+1)\top}_{i^{\prime}}\times_{\widecheck{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{k}\times_{i\in\overline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{i^{\prime}}\right)\right)_{\max(r_{k})}\right\|_{q}$
(20)

by using $\tau_{1},\ldots,\tau_{m}$, $e_{t},e_{t+1}$. The main idea to bound
(20) is to introduce
${\mathbf{I}}=P_{{\mathbf{U}}_{i^{\prime}}}+P_{{\mathbf{U}}_{i^{\prime}\perp}}$
in each mode multiplication, expand the mode products, then write the whole
term into summation of many small terms.

Step 2. After getting an upper bound for (20), we use induction to prove the
following claim,

$e_{t}\leq 2^{(d+3)/2}\tau_{1}/\lambda+e_{0}/2^{t};\quad t=0,1,\ldots.$

One technical difficulty is to deal with the sequential updating of singular
subspaces in HOOI and we use the induction idea again to tackle it. Tools we
use in this step include the singular subspace bound in [LZ20a, Theorem 6].

Step 3. The final and most challenging step involves upper bounding the tensor
reconstruction error
$\|\widetilde{{\mathbfcal{T}}}\times_{\Omega_{1}}P_{\widehat{{\mathbf{U}}}_{1}}\times\cdots\times_{\Omega_{m}}P_{\widehat{{\mathbf{U}}}_{m}}-{\mathbfcal{T}}\|_{{\rm
HS}}$ by the unified quantity $\xi$. By decomposing ${\mathbfcal{T}}$ onto the
estimated singular subspaces, we can show that

$\begin{split}&\left\|\widetilde{{\mathbfcal{T}}}\times_{\Omega_{1}}P_{\widehat{{\mathbf{U}}}_{1}}\times\cdots\times_{\Omega_{m}}P_{\widehat{{\mathbf{U}}}_{m}}-{\mathbfcal{T}}\right\|_{{\rm
HS}}\\\
\leq&\left\|{\mathbfcal{Z}}\times_{\Omega_{1}}P_{\widehat{{\mathbf{U}}}_{1}}\times\cdots\times_{\Omega_{m}}P_{\widehat{{\mathbf{U}}}_{m}}\right\|_{{\rm
HS}}+\sum_{k=1}^{d}\left\|\widehat{{\mathbf{U}}}_{k^{\prime}\perp}^{\top}\mathcal{M}_{k}({\mathbfcal{T}})\right\|_{F}.\end{split}$

By definition,
$\left\|{\mathbfcal{Z}}\times_{\Omega_{1}}P_{\widehat{{\mathbf{U}}}_{1}}\times\cdots\times_{\Omega_{m}}P_{\widehat{{\mathbf{U}}}_{m}}\right\|_{{\rm
HS}}\leq\xi$. We further show

$\begin{split}\left\|\widehat{{\mathbf{U}}}_{k^{\prime}\perp}^{\top}\mathcal{M}_{k}({\mathbfcal{T}})\right\|_{F}\overset{(a)}{\leq}C(\tau_{1},d,\lambda,e_{0})\,\xi.\end{split}$

Here $C(\tau_{1},d,\lambda,e_{0})$ is a quantity that depends on
$\tau_{1},d,\lambda,e_{0}$. The main challenge to prove (a) is that
$\widehat{{\mathbf{U}}}_{k^{\prime}}$ is not the left singular subspace of
$\mathcal{M}_{k}(\widetilde{{\mathbfcal{T}}})$. So to leverage the SVD
property of $\widehat{{\mathbf{U}}}_{k^{\prime}\perp}$, we have to project
$\widetilde{{\mathbfcal{T}}}$ onto $\widetilde{{\mathbf{U}}}_{i}^{(t_{\max})}$
and $\widetilde{{\mathbf{U}}}_{i}^{(t_{\max}-1)}$, then use the subspace
perturbation bounds established before.

Note that Theorem 3 covers the general situation where ${\mathbfcal{T}}$ may
have partial symmetric modes. We provide a corollary for the common asymmetric
case, i.e., $\Omega_{i}=\\{i\\},i=1,\ldots,d$ in the Appendix.

##  5 Perturbation Bounds of Power Iteration for Tensors with Partial Low
Multilinear Rank Structure

In some applications, e.g., multilayer network analysis [LCL19], the tensor
$\widetilde{{\mathbfcal{T}}}$ only has low-rank structure on a subset of
modes. Both the tensor power iteration algorithm and our perturbation theory
can be generalized to such cases. For a better illustration, we present the
modified algorithm and theory for order-3 tensor perturbation with mode $1$ to
be dense. Specifically, we consider

$\widetilde{{\mathbfcal{T}}}={\mathbfcal{T}}+{\mathbfcal{Z}}\in\mathbb{R}^{p_{1}\times
p_{2}\times p_{3}},$ (21)

where ${\mathbfcal{T}}$ is the signal tensor and ${\mathbfcal{Z}}$ is the
noise. We assume ${\mathbfcal{T}}$ is low-rank on modes $2$ and $3$, i.e.,
${\mathbfcal{T}}={\mathbfcal{S}}\times_{2}{\mathbf{U}}_{2}\times_{3}{\mathbf{U}}_{3}$,
where ${\mathbfcal{S}}\in\mathbb{R}^{p_{1}\times r_{2}\times r_{3}}$ is the
core tensor and ${\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}}$ for $i=2,3$. In
this setting, we consider the modified tensor power iteration algorithm for
low-rank tensor decomposition in Algorithm 3.

Algorithm 3 Power Iteration for Tensor Decomposition in Partial Multilinear
Low-Rank Setting (21)

Input: $\widetilde{{\mathbfcal{T}}}\in\mathbb{R}^{p_{1}\times p_{2}\times
p_{3}}$, initialization
$\\{\widetilde{{\mathbf{U}}}_{i}^{(0)}\in\mathbb{O}_{p_{i},r_{i}}\\}_{i=2}^{3}$,
maximum number of iterations $t_{\max}$.
Output:
$\\{\widehat{{\mathbf{U}}}_{i}\\}_{i=2}^{3},\widehat{{\mathbfcal{T}}}$.

1:For $t=1,\ldots,t_{\max}$, do
$\begin{split}&\widetilde{{\mathbf{U}}}_{2}^{(t+1)}={\rm
SVD}_{r_{2}}\left(\mathcal{M}_{2}\left(\widetilde{{\mathbfcal{T}}}\times_{3}(\widetilde{{\mathbf{U}}}_{3}^{(t)})^{\top}\right)\right)\\\
&\widetilde{{\mathbf{U}}}_{3}^{(t+1)}={\rm
SVD}_{r_{3}}\left(\mathcal{M}_{3}\left(\widetilde{{\mathbfcal{T}}}\times_{2}(\widetilde{{\mathbf{U}}}_{2}^{(t+1)})^{\top}\right)\right).\end{split}$

2:Let $\widehat{{\mathbf{U}}}_{i}=\widetilde{{\mathbf{U}}}_{i}^{(t_{\max})}$
for $i=2,3$ and compute
$\widehat{{\mathbfcal{T}}}=\widetilde{{\mathbfcal{T}}}\times_{2}P_{\widehat{{\mathbf{U}}}_{2}}\times_{3}P_{\widehat{{\mathbf{U}}}_{3}}.$

In the setting (21), we can define the blockwise errors of ${\mathbfcal{Z}}$
as follows:

$\begin{split}&\tau_{1}=\max_{k=2,3}\tau_{1k},\quad\tau_{12}=\left\|\left(\mathcal{M}_{2}({\mathbfcal{Z}}\times_{3}{\mathbf{U}}_{3}^{\top})\right)_{\max(r_{2})}\right\|_{q},\tau_{13}=\left\|\left(\mathcal{M}_{3}({\mathbfcal{Z}}\times_{2}{\mathbf{U}}_{2}^{\top})\right)_{\max(r_{3})}\right\|_{q}\\\
&\tau_{2}=\max\Big{\\{}\max_{\begin{subarray}{c}{\mathbf{V}}\in\mathbb{R}^{(p_{3}-r_{3})\times
r_{3}}\\\ \|{\mathbf{V}}\|_{q}\leq
1\end{subarray}}\left\|\left(\mathcal{M}_{2}({\mathbfcal{Z}}\times_{3}({\mathbf{U}}_{3\perp}{\mathbf{V}})^{\top})\right)_{\max(r_{2})}\right\|_{q},\\\
&\quad\quad\quad\quad\max_{\begin{subarray}{c}{\mathbf{V}}\in\mathbb{R}^{(p_{2}-r_{2})\times
r_{2}}\\\ \|{\mathbf{V}}\|_{q}\leq
1\end{subarray}}\left\|\left(\mathcal{M}_{3}({\mathbfcal{Z}}\times_{2}({\mathbf{U}}_{2\perp}{\mathbf{V}})^{\top})\right)_{\max(r_{3})}\right\|_{q}\Big{\\}};\\\
&\xi=\sup_{{\rm
rank}({\mathbfcal{Y}})\leq(p_{1},r_{2},r_{3}),\|{\mathbfcal{Y}}\|_{{\rm
HS}}\leq
1}\left\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\right\rangle.\end{split}$ (22)

We have the following perturbation bounds for the outputs of Algorithm 3.

######  Theorem 4 (Tensor Perturbation Bounds for Tensor Power Iteration with
Partial Multilinear Low-Rank (d = 3) )

Consider the perturbation model (21) with
$\widetilde{{\mathbfcal{T}}},{\mathbfcal{T}},{\mathbfcal{Z}}\in\mathbb{R}^{p_{1}\times
p_{2}\times p_{3}}$. Define the blockwise errors as in (22) and denote the
initialization error of $\\{\widetilde{{\mathbf{U}}}_{k}^{(0)}\\}_{k=2}^{3}$
as
$e_{0}:=\max_{k=2,3}\|\widetilde{{\mathbf{U}}}^{(0)\top}_{k\perp}{\mathbf{U}}_{k}\|_{q}$.
Assume the initialization error and the signal strength satisfy

$e_{0}\leq\sqrt{2}/2\,\text{ and }\,\lambda\geq 16\xi.$ (23)

Then with inputs $\widetilde{{\mathbfcal{T}}}$,
$\\{\widetilde{{\mathbf{U}}}_{i}^{(0)}\\}_{i=2}^{3}$, the mode-$k$ singular
subspace updates in Algorithm 3 after $t$ iterations satisfy

$\max_{k=2,3}\left\|\sin\Theta\left(\widetilde{{\mathbf{U}}}^{(t)}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}\leq\frac{4\sqrt{2}\tau_{1}}{\lambda}+\frac{e_{0}}{2^{t}}.$
(24)

Moreover, when $t_{\max}\geq C\log(e_{0}\lambda/\tau_{1})\vee 1$ for some
$C>0$, the outputs of estimated mode-$k$ singular subspace of Algorithm 3
satisfy

$\max_{k=2,3}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}\leq\frac{(4\sqrt{2}+1)\tau_{1}}{\lambda},$
$\begin{split}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}=\left\|\widehat{{\mathbf{U}}}_{k\perp}^{\top}{\mathbf{U}}_{k}\right\|_{q}\leq
2\sqrt{2}\left(\frac{\tau_{1k}}{\lambda}+\frac{(4\sqrt{2}+1)\tau_{1}\tau_{2}}{\lambda^{2}}\right),\quad
k=2,3.\end{split}$ (25)

The output of tensor reconstruction $\widehat{{\mathbfcal{T}}}$ satisfies

$\begin{split}\left\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\right\|_{{\rm
HS}}&\leq\left\|{\mathbfcal{Z}}\times_{2}P_{\widehat{{\mathbf{U}}}_{2}}\times_{3}P_{\widehat{{\mathbf{U}}}_{3}}\right\|_{{\rm
HS}}+\sum_{k=2}^{3}\left\|\widehat{{\mathbf{U}}}_{k\perp}^{\top}\mathcal{M}_{k}({\mathbfcal{T}})\right\|_{F}\leq(4\sqrt{2}+1)\xi.\end{split}$

The proof of Theorem 4 follows the proof of Theorem 3 easily. For simplicity,
we omit it here.

##  6 Implications in Statistics and Machine Learning

In this section, we consider a couple of applications of the HOOI perturbation
bounds we developed in statistics and machine learning. Specifically, here we
consider the perturbation model (1) and assume
${\mathbfcal{Z}}_{i_{1},\ldots,i_{d}}$’s are independent, mean-zero
$\sigma$-subgaussian, where $\sigma>0$ is the subgaussianity parameter. More
precisely,

$\mathbb{E}\exp(\lambda{\mathbfcal{Z}}_{i_{1},\ldots,i_{d}})\leq\exp(C\lambda^{2}\sigma^{2}),\,\text{for
all}\,(i_{1},\ldots,i_{d})\in[p_{1}]\times\cdots\times[p_{d}]\text{ and all
}\lambda\in\mathbb{R},$

where $C>0$ is some absolute constant. For convenience, we let
$p_{\max}=\max_{i}p_{i},p_{\min}=\min_{i}p_{i},r_{\max}=\max_{i}r_{i},r_{\min}=\min_{i}r_{i}$.

In this setting, the quantity $\xi$ is in fact closely related to the Gaussian
width [Gor88] studied in literature that measures the size or complexity of a
given set. Recall the Gaussian width of a set
$\mathcal{S}\subset\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$ is defined to
be

$w(\mathcal{S}):=\mathbb{E}\left(\sup_{{\mathbfcal{Y}}\in\mathcal{S}}\left\langle{\mathbfcal{B}},{\mathbfcal{Y}}\right\rangle\right),$

where ${\mathbfcal{B}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$ is a
tensor whose entries are independent $N(0,1)$ random variables. In view of
$w(\mathcal{S})$, we can regard $\xi$ as the Gaussian width with no
expectation and $\mathcal{S}=\\{{\mathbfcal{Y}}:\|{\mathbfcal{Y}}\|_{{\rm
HS}}\leq 1,{\rm rank}({\mathbfcal{Y}})\leq(r_{1},\ldots,r_{d})\\}$. It can be
shown in the case that ${\mathbfcal{Z}}$ has i.i.d. $N(0,1)$ entries, $\xi$
and $w(\mathcal{S})$ are the same up to constant with high probability
[RYC19].

In the following subsections, we consider two particular structures of
${\mathbfcal{T}}$, one is pure low multilinear rank structure, namely tensor
denoising/tensor PCA/tensor SVD studied in literature [RM14, ZX18, PWB16,
HSS15, LML+17, WEAM19] and another one is tensor co-clustering/block structure
[CGS+18, WZ19, HM19].

###  6.1 HOOI for Tensor Denoising

In tensor denoising, we assume ${\mathbfcal{T}}$ has the following structure,

${\mathbfcal{T}}={\mathbfcal{S}}\times_{1}{\mathbf{U}}_{1}\times\cdots\times_{d}{\mathbf{U}}_{d},$
(26)

where ${\mathbfcal{S}}\in\mathbb{R}^{r_{1}\times\cdots\times r_{d}}$ is the
core tensor and $\\{{\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}}\\}_{i=1}^{d}$
are loading matrices. With the established tensor perturbation bounds, we can
establish the following theoretical guarantee for the performance of HOOI on
tensor denoising with a very short proof.

######  Theorem 5 (Tensor Denoising: General Order $d$)

Consider the tensor denoising problem $``\eqref{eq: perturbation
model}+\eqref{eq: tensor denoising T}"$ and Algorithm 1 with inputs
$\widetilde{{\mathbfcal{T}}}$, $\Omega_{i}=\\{i\\}$, initialization
$\\{\widetilde{{\mathbf{U}}}_{i}^{(0)}\\}_{i=1}^{d}$ and
$t_{\max}=C\left(\log\frac{\lambda/\sigma}{\sqrt{p_{\max}}}\vee 1\right)$ for
some $C>0$, where
$\lambda=\min_{k}\sigma_{r_{k}}\left(\mathcal{M}_{k}({\mathbfcal{S}})\right)$
is the minimal singular value of each matricization of ${\mathbfcal{S}}$.
Assume $r_{\max}\leq p_{\min}^{1/(d-1)}$ and
$\max_{i}\|\widetilde{{\mathbf{U}}}_{i\perp}^{(0)\top}{\mathbf{U}}_{i}\|\leq\sqrt{2}/2$.
Then if $\lambda/\sigma\geq
2^{(d+4)/2}(1+\sqrt{2}/2)^{d-1}\sqrt{p_{\max}r_{\max}}$, with probability at
least $1-\exp(-cp_{\min})$, the output
$\widehat{{\mathbf{U}}}_{k},\widehat{{\mathbfcal{T}}}$ satisfy

$\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|\leq
C\left(1-\frac{c}{r_{\max}}\right)^{-(d-1)/2}\left(\frac{\sqrt{p_{k}}}{\lambda/\sigma}+\frac{\sqrt{p_{\max}^{2}r_{\max}}}{\left(\lambda/\sigma\right)^{2}}\right)$

and

$\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm HS}}\leq
C\left(1+2d\left(1-\frac{c}{r_{\max}}\right)^{-\frac{d-1}{2}}\right)\sigma\sqrt{\sum_{i=1}^{d}p_{i}r_{i}},$
(27)

for some constants $c,C>0$.

When $d$ is a constant, the upper bound for tensor reconstruction error
matches the lower bound in [ZX18, Theorem 3], which shows HOOI achieves the
optimal tensor reconstruction error in the tensor denoising problem.

###  6.2 Tensor Co-clustering/Block Model

Co-clustering is among the most important unsupervised learning methods that
reveals the checkerbox-like association pattern in data. A number of
algorithms have been proposed [WLP15, WBG16, KS08, PSB12, ZWC+16, SPL+09,
JSB09] for tensor co-clustering in the literature, however most of the work
does not provide statistical guarantees for recovering the underlying co-
clustering structure. Very recently, [CGS+18] and [WZ19] studied the
performance of co-clustering estimation and cocluster recovery based on convex
relaxation and combinatorial search algorithms. By using the tools of
perturbation bounds of HOOI given in Section 4, we are able to provide the
first guarantee for co-clustering estimation and cocluster recovery based on
computational efficient HOOI algorithm. Compared to the convex relaxation
approach [CGS+18], HOOI has a better guarantee for tensor reconstruction and
it also gives guarantee for cocluster membership recovery. Specifically, in
the tensor co-clustering/block model, we assume ${\mathbfcal{T}}$ has the
following structure,

${\mathbfcal{T}}={\mathbfcal{B}}\times_{1}{\mathbf{\Pi}}_{1}\times\cdots\times_{d}{\mathbf{\Pi}}_{d},$
(28)

where ${\mathbf{\Pi}}_{i}\in\mathbb{M}_{p_{i},r_{i}}$ and
$\mathbb{M}_{p_{i},r_{i}}$ is the collection of all $p_{i}\times r_{i}$
membership matrices with each row has exactly one $1$ and $(r_{i}-1)$ $0$’s.
For any ${\mathbf{\Pi}}_{i}$, the cocluster membership of node $k$ is denoted
by $g_{i}^{(k)}\in[r_{i}]$, which satisfies
$({\mathbf{\Pi}}_{i})_{[k,g_{i}^{(k)}]}=1$. Let $G^{(j)}_{i}\equiv
G^{(j)}_{i}({\mathbf{\Pi}}_{i})=\\{k\in[p_{i}]:g_{i}^{(k)}=j\\}$ be the set of
$i$th mode node indices that belongs to cocluster $j$ and
$p^{(j)}_{i}=|G^{(j)}_{i}|$ for all $j\in[r_{i}]$. For simplicity, we assume
the cocluster sizes for each cluster are on the same order for every mode,
i.e.,

$p_{i}^{(1)}\asymp p_{i}^{(2)}\asymp\cdots\asymp
p_{i}^{(r_{i})}\asymp\frac{p_{i}}{r_{i}},\text{ for }i=1,\ldots,d.$ (A1)

We consider two cocluster membership recovery error metrics:

  1. 1.

Let $\|{\mathbf{M}}\|_{0}$ be the number of nonzero entries in matrix
${\mathbf{M}}$. Suppose ${\mathbf{E}}_{r_{i}}$ is the set of all $r_{i}\times
r_{i}$ permutation matrices. Define

$err(\widehat{{\mathbf{\Pi}}}_{i},{\mathbf{\Pi}}_{i})=\frac{1}{p_{i}}\min_{{\mathbf{J}}\in{\mathbf{E}}_{r_{i}}}\|\widehat{{\mathbf{\Pi}}}_{i}{\mathbf{J}}-{\mathbf{\Pi}}_{i}\|_{0}$
(29)

as the misclassification rate of $\widehat{{\mathbf{\Pi}}}_{i}$.

  2. 2.

Define

$\widetilde{err}(\widehat{{\mathbf{\Pi}}}_{i},{\mathbf{\Pi}}_{i})=\min_{{\mathbf{J}}\in{\mathbf{E}}_{r_{i}}}\max_{1\leq
j\leq
r_{i}}\frac{1}{p^{(j)}_{i}}\|(\widehat{{\mathbf{\Pi}}}_{i}{\mathbf{J}})_{[G^{(j)}_{i},:]}-({\mathbf{\Pi}}_{i})_{[G^{(j)}_{i},:]}\|_{0}.$

Intuitively speaking,
$\widetilde{err}(\widehat{{\mathbf{\Pi}}}_{i},{\mathbf{\Pi}}_{i})$ measures
the worst relative misclassification rates over all communities.

It is easy to check that $0\leq
err(\widehat{{\mathbf{\Pi}}}_{i},{\mathbf{\Pi}}_{i})\leq\widetilde{err}(\widehat{{\mathbf{\Pi}}}_{i},{\mathbf{\Pi}}_{i})\leq
2$.

The following lemma gives a Tucker decomposition of ${\mathbfcal{T}}$ in the
tensor block model (28). This decomposition bridges the tensor co-clustering
model and the tensor signal-plus-noise model, which explains why HOOI would
work for tensor co-clustering.

###### Lemma 3

Suppose ${\mathbfcal{T}}$ has the tensor co-clustering/block structure (28),
where ${\mathbfcal{B}}$ is a multilinear rank-$(r_{1},\ldots,r_{d})$ tensor.
Assume
${\mathbfcal{S}}\times_{1}{\mathbf{V}}_{1}\times\cdots\times_{d}{\mathbf{V}}_{d}$
with ${\mathbf{V}}_{i}\in\mathbb{O}_{r_{i},r_{i}}$ is the Tucker decomposition
of
${\mathbfcal{B}}\times_{1}({\mathbf{\Pi}}_{1}^{\top}{\mathbf{\Pi}}_{1})^{\frac{1}{2}}\times\cdots\times({\mathbf{\Pi}}_{d}^{\top}{\mathbf{\Pi}}_{d})^{\frac{1}{2}}$.
Then

${\mathbfcal{T}}={\mathbfcal{S}}\times_{1}{\mathbf{U}}_{1}\times\cdots\times_{d}{\mathbf{U}}_{d}$

with
${\mathbf{U}}_{i}={\mathbf{\Pi}}_{i}({\mathbf{\Pi}}_{i}^{\top}{\mathbf{\Pi}}_{i})^{-\frac{1}{2}}{\mathbf{V}}_{i}\in\mathbb{O}_{p_{i},r_{i}}$
for $i=1,\ldots,d$.

The following Theorem 6 gives the theoretical guarantee on the performance of
Algorithm 4 for tensor reconstruction and cocluster membership recovery.

Algorithm 4 HOOI for Tensor Co-clustering/Block Model

Input: Tensor
$\widetilde{{\mathbfcal{T}}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$,
indices group $\Omega_{i}=\\{i\\}$, initialization
$\widetilde{{\mathbf{U}}}^{(0)}_{i}\in\mathbb{R}^{p_{i}\times r_{i}}$ for
$i=1,\ldots,d$, maximum number of iterations $t_{\max}$,
Output:
$\widehat{{\mathbf{\Pi}}}_{i}\in\mathbb{M}_{p_{i},r_{i}}$,$i=1,\ldots,d$ and
$\widehat{{\mathbfcal{T}}}$.

1:Apply Algorithm 1 with input $\widetilde{{\mathbfcal{T}}}$,
$\\{\Omega_{i}\\}_{i=1}^{d}$,
$\\{\widetilde{{\mathbf{U}}}^{(0)}_{i}\\}_{i=1}^{d}$, maximum number of
iterations $t_{\max}$ and get outputs
$\\{\widehat{{\mathbf{U}}}_{i}\\}_{i=1}^{d}$ and $\widehat{{\mathbfcal{T}}}$.

2:For each mode $i$, apply $\epsilon$-approximation K-means [KSS04] on
$\widehat{{\mathbf{U}}}_{i}$, i.e., compute
$\widehat{{\mathbf{\Pi}}}_{i},\widehat{{\mathbf{X}}}_{i}\in\mathbb{M}_{p_{i},r_{i}}\times\mathbb{R}^{r_{i}\times
r_{i}}$ such that
$\|\widehat{{\mathbf{\Pi}}}_{i}\widehat{{\mathbf{X}}}_{i}-\widehat{{\mathbf{U}}}_{i}\|_{F}^{2}\leq(1+\epsilon)\min_{{\mathbf{\Pi}},{\mathbf{X}}\in\mathbb{M}_{p_{i},r_{i}}\times\mathbb{R}^{r_{i}\times
r_{i}}}\|{\mathbf{\Pi}}{\mathbf{X}}-\widehat{{\mathbf{U}}}_{i}\|_{F}^{2}.$
(30)

######  Theorem 6 (HOOI for Tensor Co-clustering/Block Model)

Consider the tensor co-clustering/block model “$\eqref{eq: perturbation
model}+\eqref{eq: tensor block model}$” and the Algorithm 4 with inputs
$\widetilde{{\mathbfcal{T}}}$, initializations
$\\{\widetilde{{\mathbf{U}}}_{i}^{(0)}\\}_{i=1}^{d}$ and
$t_{\max}=C\left(\log\left(\frac{\lambda/\sigma}{\sqrt{p_{\max}}}\sqrt{\frac{\prod_{i=1}^{d}p_{i}}{\prod_{i=1}^{d}r_{i}}}\right)\vee
1\right)$ for some $C>0$, where
$\lambda=\min_{k}\sigma_{r_{k}}\left(\mathcal{M}_{k}({\mathbfcal{B}})\right)$
is the minimal singular value at each matricization of the core tensor
parameter ${\mathbfcal{B}}$. Assume $r_{\max}\leq p_{\min}^{1/(d-1)}$,
$\max_{i}\|\widetilde{{\mathbf{U}}}^{(0)\top}_{i\perp}{\mathbf{U}}_{i}\|\leq\frac{\sqrt{2}}{2}$,
and (A1) holds. Then if

$\lambda/\sigma\geq
C\sqrt{\frac{p_{\max}^{2}r_{\max}\prod_{i=1}^{d}r_{i}}{p_{\min}\prod_{i=1}^{d}p_{i}}},$

for sufficiently large constant $C>0$, with probability at least
$1-\exp(-cp_{\min})$ for some $c>0$, $\widehat{{\mathbf{U}}}_{k}$,
$\widehat{{\mathbfcal{T}}}$ satisfy

$\|\sin\Theta(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k})\|\leq
C(d)\frac{\sqrt{p_{k}}}{\lambda/\sigma}\sqrt{\frac{\prod_{i=1}^{d}r_{i}}{\prod_{i=1}^{d}p_{i}}},\quad
k=1,\ldots,d$ $\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm HS}}\leq
C(d)\sigma\sqrt{\sum_{i=1}^{d}p_{i}r_{i}};$

and we also have the following upper bound on cocluster recovery error,

$\begin{split}&err(\widehat{{\mathbf{\Pi}}}_{i},{\mathbf{\Pi}}_{i})\leq
C_{1}(d,\epsilon)\frac{p^{\prime}_{i,\max}r_{i}}{(\lambda/\sigma)^{2}}\frac{\prod_{i=1}^{d}r_{i}}{\prod_{i=1}^{d}p_{i}},\\\
&\widetilde{err}(\widehat{{\mathbf{\Pi}}}_{i},{\mathbf{\Pi}}_{i})\leq
C_{2}(d,\epsilon)\frac{p_{i}r_{i}}{(\lambda/\sigma)^{2}}\frac{\prod_{i=1}^{d}r_{i}}{\prod_{i=1}^{d}p_{i}}.\end{split}$

Here $C(d),C_{1}(d,\epsilon),C_{2}(d,\epsilon)>0$ are some constants depending
only on $d$ and $\epsilon$, $p^{\prime}_{i,\max}$ is the second largest
cocluster size at mode $i$.

Note that our cocluster recovery guarantee is new for polynomial-time
algorithms. When $p=p_{1}=\cdots=p_{d}$, the best polynomial time algorithm
guarantee for tensor reconstruction is $p^{d-1}$ in [CGS+18] and our result
can be significantly better.

##  7 Numerical Studies

In this section, we first provide numerical studies to support the main
theoretical results in Section 4 and then compare HOOI with other existing
algorithms for tensor decomposition in applications of tensor denoising and
tensor co-clustering. Throughout the simulation, we consider order-$3$ tensor
perturbation setting
$\widetilde{{\mathbfcal{T}}}={\mathbfcal{T}}+{\mathbfcal{Z}}$ with
${\mathbfcal{Z}}$ being the noise tensor with i.i.d. $N(0,\sigma^{2})$
entries. Without particular specification, we set
$p=p_{1}=p_{2}=p_{3},r=r_{1}=r_{2}=r_{3}$. The error metrics we consider for
tensor reconstruction and mode-$k$ singular subspace estimation are root mean
square error (RMSE) $\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm HS}}$
and $\|\sin\Theta(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k})\|$,
respectively. All simulations are repeated $100$ times and the average
statistics are reported.

###  7.1 Perturbation Bounds of HOOI with good initialization

In this simulation, we study the perturbation bounds of HOOI with randomly
generated good initialization. Let
${\mathbfcal{T}}={\mathbfcal{S}}\times_{1}{\mathbf{U}}_{1}\times_{2}{\mathbf{U}}_{2}\times_{3}{\mathbf{U}}_{3}$,
where ${\mathbf{U}}_{i}\in\mathbb{R}^{p_{i}\times r}$ is generated uniformly
at random from $\mathbb{O}_{p_{i},r}$ and
${\mathbfcal{S}}\in\mathbb{R}^{r\times r\times r}$ is a diagonal tensor with
diagonal values $\\{i\lambda\\}_{i=1}^{r}$. The initializations of
${\mathbf{U}}_{i}$ of Algorithm 1 are
$\widetilde{{\mathbf{U}}}_{i}^{(0)}=\frac{1}{\sqrt{2}}{\mathbf{U}}_{i}+\frac{1}{\sqrt{2}}{\mathbf{U}}_{i}^{\prime}$,
where ${\mathbf{U}}_{i}^{\prime}={\mathbf{U}}_{i^{\prime}\perp}\mathbf{O}$ for
some random orthogonal matrix $\mathbf{O}\in\mathbb{O}_{p_{i}-r,r}$. It is
easy to check that
$\|\sin\Theta({\mathbf{U}}_{i},\widetilde{{\mathbf{U}}}_{i}^{(0)})\|=\frac{\sqrt{2}}{2}$
for $i=1,2,3$.

First for tensor reconstruction, let $p\in\\{20,30,\ldots,100\\},r=5$,
$\sigma\in\\{1,2,3,4\\}$ and $\lambda=5\sqrt{pr}\sigma$. We can check that
with high probability, $\|{\mathbfcal{Z}}\|_{{\rm HS}}\leq
Cp^{\frac{3}{2}}\sigma$ and $\xi\leq C\sqrt{pr}\sigma$ for some $C>0$
following the same proof as Theorem 5. In Figure 2(a), the RMSE of tensor
reconstruction of HOOI is presented. We find as the perturbation results in
Section 4 suggest, $\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm HS}}$
can be much smaller than $\|{\mathbfcal{Z}}\|_{{\rm HS}}$. This demonstrates
the superior performance of the HOOI estimator compared to the trivial
estimator $\widetilde{{\mathbfcal{T}}}$. At the same time, the RMSE for tensor
reconstruction increases as $p$ and $\sigma$ become bigger and this matches
our theoretical findings in Theorem 3 that the error bound of HOOI for
$\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm HS}}$ is $O(\xi)$, which
increases as $p,\sigma$ increase.

(a) Tensor Reconstruction (b) Convergence of Singular Subspace

Figure 2: HOOI with good initialization. (a) Tensor reconstruction error
$\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm HS}}$ for
$p\in\\{20,30,\ldots,100\\},r=5$, $\sigma\in\\{1,2,3,4\\}$ and
$\lambda=5\sqrt{pr}\sigma$; (b) Mode-$k$ singular subspace estimation with and
without rescaling under $p_{1}=10,p_{2}=100,p_{3}=500$,
$r\in\\{3,5\\},\sigma=1$ and $\lambda=\alpha\cdot
p_{3}\frac{\sqrt{r_{1}}}{\sqrt{p_{1}}}$ with varying $\alpha$

Next we demonstrate the unilateral perturbation bounds for mode-$k$ singular
subspace estimation. Specifically, we consider $p_{1}=10,p_{2}=100,p_{3}=500$,
$r\in\\{3,5\\},\sigma=1$ and $\lambda=\alpha\cdot
p_{3}\frac{\sqrt{r}}{\sqrt{p_{1}}}$ with varying $\alpha$. The errors of the
mode-1, mode-2, mode-3 estimated singular subspaces with and without rescaling
are provided in Figure 2(b). We can see from Figure 2(b) left panel the errors
of estimated singular subspaces converge to different values depending on the
corresponding mode size $p_{i}$, and a further rescaling of estimation error
by $\sqrt{p_{i}}$ makes them roughly on the same level (see Figure 2(b) right
panel). This matches the unilateral property of the singular subspace
perturbation results in Remark 3 that when
$\lambda=O(p_{3}\frac{\sqrt{r}}{\sqrt{p_{1}}})$,
$\|\sin\Theta(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k})\|\leq
C\frac{\sqrt{p_{i}}}{\lambda/\sigma}$ for some $C>0$, and this upper bound
increases linearly with respect to $\sqrt{p_{i}}$.

###  7.2 Comparison of HOOI with other Algorithms

In this section, we do a comparison of HOOI with truncated HOSVD (T-HOSVD)
[DLDMV00a] and sequentially truncated HOSVD (ST-HOSVD) [VVM12] in the tensor
denoising and tensor co-clustering applications. We also include one-step HOOI
(O-HOOI), since it might be useful as a surrogate of HOOI in large scale
tensor decomposition settings as we mentioned in Remark 5. The initialization
we consider for HOOI and O-HOOI are ST-HOSVD with natural truncation order,
i.e., $\widetilde{{\mathbf{U}}}^{(0)}_{i}={\rm
SVD}_{r}(\mathcal{M}_{i}(\widetilde{{\mathbfcal{T}}}\times_{j<i}\widetilde{{\mathbf{U}}}_{j}^{(0)}))$.
In Table 1, we give the time complexity of HOOI, O-HOOI, ST-HOSVD and T-HOSVD.
We can see that as long as $dr\leq p$, a common case in practice, the time
complexity of O-HOOI and ST-HOSVD are on the same order, and they could be
faster than full HOOI and HOSVD in general.

| HOOI | O-HOOI | ST-HOSVD | T-HOSVD
---|---|---|---|---
Complexity | $O(p^{d+1}+t_{\max}drp^{d})$ | $O(p^{d+1}+drp^{d})$ | $O(p^{d+1})$ | $O(dp^{d+1})$
Table 1: Time Complexity of HOOI, O-HOOI, ST-HOSVD and T-HOSVD under setting
$p_{1}=\cdots=p_{d}=p,r_{1}=\cdots=r_{d}=r$, $r\ll p$. HOOI and O-HOOI is
initialized by ST-HOSVD.

In tensor denoising, the generating process of ${\mathbfcal{T}}$ is the same
as before. Let $p=100$, $r=5$, $\sigma\in\\{1,2\\}$, $\lambda=\alpha\cdot
p^{\frac{3}{4}}\sigma$ with varying $\alpha$. The comparison of these
algorithms for tensor reconstruction and singular subspace estimation are
given in Figure 3. First, we find that HOOI is best in both tensor
reconstruction and singular subspace estimation among four algorithms.
Meanwhile, O-HOOI is slightly worse than HOOI for small $\alpha$ and has very
close performance with HOOI when $\alpha$ is relative large, which suggests
that in some computationally heavy applications, we can just run HOOI for one
iteration to achieve reasonable estimation. Part of this phenomenon can be
explained by the one-iteration optimality of HOOI for tensor reconstruction as
we discussed in Remark 5. On the other hand, HOOI and O-HOOI are often much
better than T-HOSVD and ST-HOSVD for both tensor reconstruction and mode-$k$
singular subspace estimation within a wide range of settings.

(a) Tensor Reconstruction (b) Singular Subspace Estimation

Figure 3: Comparison of HOOI, one-step HOOI (O-HOOI), truncated HOSVD
(T-HOSVD), sequentially truncated HOSVD (ST-HOSVD) in tensor denoising under
$p=100$, $r=5$, $\sigma\in\\{1,2\\}$, $\lambda=\alpha\cdot
p^{\frac{3}{4}}\sigma$ with $\alpha\in[1,4]$. (a) Tensor reconstruction; (b)
Averaged singular subspace estimation

Finally, we study the performance of the HOOI-based Algorithm 4 in tensor co-
clustering recovery and do a comparison of it with T-HOSVD, ST-HOSVD, and
O-HOOI based clustering algorithms. In this simulation, we generate
${\mathbfcal{T}}={\mathbfcal{B}}\times_{1}{\mathbf{\Pi}}_{1}\times_{2}{\mathbf{\Pi}}_{2}\times_{3}{\mathbf{\Pi}}_{3}$
such that $\\{{\mathbf{\Pi}}_{i}\\}_{i=1}^{3}$ have balanced cluster size and
${\mathbfcal{B}}=\frac{{\mathbfcal{B}}_{0}}{\min_{i}\sigma_{r}\left(\mathcal{M}_{i}({\mathbfcal{B}}_{0})\right)}\lambda$
with ${\mathbfcal{B}}_{0}\overset{i.i.d.}{\sim}N(0,1)$. The error metric we
consider is the average cocluster membership misclassification error rate in
(29). The performance of Algorithm 4 under
$p\in\\{50,80\\},r\in\\{3,5,8\\},\sigma=1,\lambda=\alpha\cdot\frac{r^{3/2}}{p^{3/4}}\sigma$
is presented in Figure 4(a). We can see the misclassification error decreases
as the signal strength increases and cocluster number decreases. The
comparison of Algorithm 4 and T-HOSVD, ST-HOSVD, O-HOOI based spectral
clustering is given in Figure 4(b) under the same setting with $r=5$. Again,
HOOI-based algorithm has the best performance in cocluster recovery. O-HOOI
and ST-HOSVD perform similarly here and both of them are much better than
T-HOSVD.

(a) HOOI for Cocluster Recovery (b) Comparison in Cocluster Recovery

Figure 4: Tensor cocluster recovery under $\sigma=1$,
$\lambda=\alpha\cdot\frac{r^{3/2}}{p^{3/4}}\sigma$ with varying $\alpha$. (a)
HOOI on tensor cocluster recovery under $p\in\\{50,80\\}$, $r\in\\{3,5,8\\}$.
(b) Comparison of HOOI, one-step HOOI (O-HOOI), truncated HOSVD (T-HOSVD),
sequentially truncated HOSVD (ST-HOSVD) in cocluster recovery under
$p=80,r=5$.

##  8 Conclusion and Discussion

In this paper, we provide the first sharp blockwise perturbation bounds of
HOOI for tensors with guarantees for both tensor reconstruction and mode-$k$
singular subspace estimation. Furthermore, we show both HOOI and one-step HOOI
with good initialization is optimal in terms of tensor reconstruction by
providing rate matching lower bound. Finally, we support our theoretical
results with extensive numerical studies and apply them to tensor denoising
and tensor co-clustering applications. Apart from the applications mentioned
above, the main perturbation results can be applied to many other applications
where tensor “spectral method” HOOI is applicable, such as tensor completion
[YZ14, YZ17, XYZ17, XY17], hypergraphic stochastic block model [GD14, GD17,
KSX19, CLW19, ALS18, KBG18], multilayer network [LCL19, JLLX20], MPCA [LPV08],
latent variable model [AGH+14], etc. In tensor completion [XYZ17] and many
other applications, more specialized initializers can achieve better
performance than HOSVD – the classic initializers for HOOI in the literature
[DLDMV00b]. Our tensor perturbation bounds still apply to these cases as our
theoretical analysis admits all initializers satisfying certain mild
conditions.

At the same time, due to the NP hardness of computing many tensor quantities
[HL13], the Alternating Least Square (ALS) and Power iteration have been the
“workhorse” algorithms in computing low-rank tensor approximation and solving
many other tensor problems [KB09]. Our induction proof idea in Theorem 3 could
also shed light on how to analyze other iterative ALS/Power iteration
procedures for tensor problems [ZLZ13, LPV08, WZ19, XYZ+05, YXY+05, YXY+06,
LXZL18, LW20].

Finally, an interesting open question is to provide perturbation results for
other algorithms, such as (quasi-)Newton-Grassmann method [ES09, SL10],
geometric Newton method [IDLAVH09], Riemannian trust region scheme [IAVHDL11]
since sometimes Riemannian trust region scheme and/or Newton-type methods can
take much fewer iterations than HOOI to converge. Also in the paper, we mainly
focus on Tucker format of tensor decomposition. Although it has many
advantages, it has the drawback that in ultra high-order tensor problems, the
storage cost of the core tensor in Tucker format scales exponentially with
respect to the tensor order and it is more desirable to consider other low-
rank tensor approximation methods than Tucker, such as Hierarchical Tucker
decomposition [BG13, Gra10, HK09], and Tensor Train decomposition [Ose11,
OT09]. It is interesting to develop the perturbation bounds for algorithms on
Hierarchical Tucker or Tensor Train tensor decomposition.

## References

  * [AFWZ17] Emmanuel Abbe, Jianqing Fan, Kaizheng Wang, and Yiqiao Zhong.  Entrywise eigenvector analysis of random matrices with low expected rank.  arXiv preprint arXiv:1709.09565, 2017.
  * [AGH+14] Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M Kakade, and Matus Telgarsky.  Tensor decompositions for learning latent variable models.  Journal of Machine Learning Research, 15:2773–2832, 2014.
  * [AGJ14] Animashree Anandkumar, Rong Ge, and Majid Janzamin.  Guaranteed non-orthogonal tensor decomposition via alternating rank-$1$ updates.  arXiv preprint arXiv:1402.5180, 2014.
  * [AJSN16] Anima Anandkumar, Prateek Jain, Yang Shi, and Uma Naresh Niranjan.  Tensor vs. matrix methods: Robust tensor decomposition under block sparse perturbations.  In Artificial Intelligence and Statistics, pages 268–276, 2016\.
  * [ALS18] Kwangjun Ahn, Kangwook Lee, and Changho Suh.  Hypergraph spectral clustering in the weighted stochastic block model.  IEEE Journal of Selected Topics in Signal Processing, 12(5):959–974, 2018.
  * [AY20] Arnab Auddy and Ming Yuan.  Perturbation bounds for orthogonally decomposable tensors and their applications in high dimensional data analysis.  arXiv preprint arXiv:2007.09024, 2020.
  * [BCS+16] Emilio Rafael Balda, Sher Ali Cheema, Jens Steinwandt, Martin Haardt, Amir Weiss, and Arie Yeredor.  First-order perturbation analysis of low-rank tensor approximations based on the truncated hosvd.  In 2016 50th Asilomar Conference on Signals, Systems and Computers, pages 1723–1727. IEEE, 2016.
  * [BFY18] Krishnakumar Balasubramanian, Jianqing Fan, and Zhuoran Yang.  Tensor methods for additive index models under discordance and heterogeneity.  arXiv preprint arXiv:1807.06693, 2018.
  * [BG13] Jonas Ballani and Lars Grasedyck.  A projection method to solve linear systems in tensor format.  Numerical linear algebra with applications, 20(1):27–43, 2013.
  * [BGN11] Florent Benaych-Georges and Raj Rao Nadakuditi.  The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices.  Advances in Mathematics, 227(1):494–521, 2011.
  * [BGN12] Florent Benaych-Georges and Raj Rao Nadakuditi.  The singular values and vectors of low rank perturbations of large rectangular random matrices.  Journal of Multivariate Analysis, 111:120–135, 2012.
  * [Bha13] Rajendra Bhatia.  Matrix analysis, volume 169.  Springer Science & Business Media, 2013.
  * [BK12] Brett W Bader and Tamara G Kolda.  Matlab tensor toolbox version 2.5.  Available online, January, 7, 2012.
  * [BQS18] Xuan Bi, Annie Qu, and Xiaotong Shen.  Multilayer tensor factorization with applications to recommender systems.  The Annals of Statistics, 46(6B):3308–3333, 2018.
  * [CC10] Cesar F Caiafa and Andrzej Cichocki.  Generalizing the column–row matrix decomposition to multi-way arrays.  Linear Algebra and its Applications, 433(3):557–573, 2010.
  * [CCF18] Yuxin Chen, Chen Cheng, and Jianqing Fan.  Asymmetry helps: Eigenvalue and eigenvector analyses of asymmetrically perturbed low-rank matrices.  arXiv preprint arXiv:1811.12804, 2018.
  * [CCT12] Kamalika Chaudhuri, Fan Chung, and Alexander Tsiatas.  Spectral clustering of graphs with general degrees in the extended planted partition model.  In Conference on Learning Theory, pages 35–1, 2012.
  * [CCZ16] Tianxi Cai, T. Tony Cai, and Anru Zhang.  Structured matrix completion with applications to genomic data integration.  Journal of the American Statistical Association, 111(514):621–633, 2016.
  * [CGS+18] 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.
  * [CLW19] I Chien, Chung-Yi Lin, and I-Hsiang Wang.  On the minimax misclassification ratio of hypergraph community detection.  IEEE Transactions on Information Theory, 2019.
  * [CMDL+15] Andrzej Cichocki, Danilo Mandic, Lieven De Lathauwer, Guoxu Zhou, Qibin Zhao, Cesar Caiafa, and Huy Anh Phan.  Tensor decompositions for signal processing applications: From two-way to multiway component analysis.  IEEE Signal Processing Magazine, 32(2):145–163, 2015.
  * [CRV15] Peter Chin, Anup Rao, and Van Vu.  Stochastic block model and community detection in sparse graphs: A spectral algorithm with optimal rate of recovery.  In Conference on Learning Theory, pages 391–423, 2015.
  * [CRY16] Han Chen, Garvesh Raskutti, and Ming Yuan.  Non-convex projected gradient descent for generalized low-rank tensor regression.  arXiv preprint arXiv:1611.10349, 2016.
  * [CTP19a] Joshua Cape, Minh Tang, and Carey E Priebe.  Signal-plus-noise matrix models: eigenvector deviations and fluctuations.  Biometrika, 106(1):243–250, 2019.
  * [CTP19b] Joshua Cape, Minh Tang, and Carey E Priebe.  The two-to-infinity norm and singular subspace geometry with applications to high-dimensional statistics.  The Annals of Statistics, 47(5):2405–2439, 2019.
  * [CZ18] T Tony Cai and Anru Zhang.  Rate-optimal perturbation bounds for singular subspaces with applications to high-dimensional statistics.  The Annals of Statistics, 46(1):60–89, 2018.
  * [DK70] Chandler Davis and William Morton Kahan.  The rotation of eigenvectors by a perturbation. iii.  SIAM Journal on Numerical Analysis, 7(1):1–46, 1970.
  * [DL04] Lieven De Lathauwer.  First-order perturbation analysis of the best rank-(r1, r2, r3) approximation in multilinear algebra.  Journal of Chemometrics: A Journal of the Chemometrics Society, 18(1):2–11, 2004.
  * [DLDMV00a] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle.  A multilinear singular value decomposition.  SIAM journal on Matrix Analysis and Applications, 21(4):1253–1278, 2000.
  * [DLDMV00b] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle.  On the best rank-1 and rank-(r 1, r 2,…, rn) approximation of higher-order tensors.  SIAM journal on Matrix Analysis and Applications, 21(4):1324–1342, 2000.
  * [DLV04] Lieven De Lathauwer and Joos Vandewalle.  Dimensionality reduction in higher-order signal processing and rank-(r1, r2,, rn) reduction in multilinear algebra.  Linear Algebra and its Applications, 391:31–55, 2004.
  * [Dop00] Froilán M Dopico.  A note on sin $\theta$ theorems for singular subspace variations.  BIT Numerical Mathematics, 40(2):395–403, 2000.
  * [EBW17] Justin Eldridge, Mikhail Belkin, and Yusu Wang.  Unperturbed: spectral analysis beyond davis-kahan.  arXiv preprint arXiv:1706.06516, 2017.
  * [ES09] Lars Elden and Berkant Savas.  A newton–grassmann method for computing the best multilinear rank-(r1, r2, r3) approximation of a tensor.  SIAM Journal on Matrix Analysis and applications, 31(2):248–271, 2009.
  * [ES11] Lars Elden and Berkant Savas.  Perturbation theory and optimality conditions for the best multilinear rank approximation of a tensor.  SIAM journal on matrix analysis and applications, 32(4):1422–1450, 2011.
  * [EY36] Carl Eckart and Gale Young.  The approximation of one matrix by another of lower rank.  Psychometrika, 1(3):211–218, 1936.
  * [FWZ16] Jianqing Fan, Weichen Wang, and Yiqiao Zhong.  An l infinity eigenvector perturbation bound and its application to robust covariance estimation.  arXiv preprint arXiv:1603.03516, 2016.
  * [GD14] Debarghya Ghoshdastidar and Ambedkar Dukkipati.  Consistency of spectral partitioning of uniform hypergraphs under planted partition model.  In Advances in Neural Information Processing Systems, pages 397–405, 2014.
  * [GD17] Debarghya Ghoshdastidar and Ambedkar Dukkipati.  Consistency of spectral hypergraph partitioning under planted partition model.  The Annals of Statistics, 45(1):289–315, 2017.
  * [GGH14] Quanquan Gu, Huan Gui, and Jiawei Han.  Robust tensor decomposition with gross corruption.  In Advances in Neural Information Processing Systems, pages 1422–1430, 2014.
  * [GHS87] Gene H Golub, Alan Hoffman, and Gilbert W Stewart.  A generalization of the eckart-young-mirsky matrix approximation theorem.  Linear Algebra and its applications, 88:317–327, 1987.
  * [GKT13] Lars Grasedyck, Daniel Kressner, and Christine Tobler.  A literature survey of low-rank tensor approximation techniques.  GAMM-Mitteilungen, 36(1):53–78, 2013.
  * [Gor88] Yehoram Gordon.  On milman’s inequality and random subspaces which escape through a mesh in ? n.  In Geometric aspects of functional analysis, pages 84–106. Springer, 1988.
  * [GQ14] Donald Goldfarb and Zhiwei Qin.  Robust low-rank tensor recovery: Models and algorithms.  SIAM Journal on Matrix Analysis and Applications, 35(1):225–253, 2014.
  * [Gra10] Lars Grasedyck.  Hierarchical singular value decomposition of tensors.  SIAM Journal on Matrix Analysis and Applications, 31(4):2029–2054, 2010.
  * [HK09] Wolfgang Hackbusch and Stefan Kühn.  A new scheme for the tensor representation.  Journal of Fourier analysis and applications, 15(5):706–722, 2009\.
  * [HL13] Christopher J Hillar and Lek-Heng Lim.  Most tensor problems are np-hard.  Journal of the ACM (JACM), 60(6):45, 2013.
  * [HM19] Rui Henriques and Sara C Madeira.  Triclustering algorithms for three-dimensional data analysis: A comprehensive survey.  ACM Computing Surveys (CSUR), 51(5):95, 2019.
  * [HSS15] Samuel B Hopkins, Jonathan Shi, and David Steurer.  Tensor principal component analysis via sum-of-square proofs.  In Conference on Learning Theory, pages 956–1006, 2015.
  * [HWZ20] Rungang Han, Rebecca Willett, and Anru Zhang.  An optimal statistical and computational framework for generalized tensor estimation.  arXiv preprint arXiv:2002.11255, 2020.
  * [HZC18] Botao Hao, Anru Zhang, and Guang Cheng.  Sparse and low-rank tensor estimation via cubic sketchings.  arXiv preprint arXiv:1801.09326, 2018.
  * [IAVHDL11] Mariya Ishteva, P-A Absil, Sabine Van Huffel, and Lieven De Lathauwer.  Best low multilinear rank approximation of higher-order tensors, based on the riemannian trust-region scheme.  SIAM Journal on Matrix Analysis and Applications, 32(1):115–135, 2011.
  * [IDLAVH09] Mariya Ishteva, Lieven De Lathauwer, P-A Absil, and Sabine Van Huffel.  Differential-geometric newton method for the best rank-(r1, r2,r3) approximation of tensors.  Numerical Algorithms, 51(2):179–194, 2009.
  * [JLLX20] Bing-Yi Jing, Ting Li, Zhongyuan Lyu, and Dong Xia.  Community detection on mixture multi-layer networks via regularized tensor decomposition.  arXiv preprint arXiv:2002.04457, 2020.
  * [JSB09] Stefanie Jegelka, Suvrit Sra, and Arindam Banerjee.  Approximation algorithms for tensor clustering.  In International Conference on Algorithmic Learning Theory, pages 368–383. Springer, 2009.
  * [KB09] Tamara G Kolda and Brett W Bader.  Tensor decompositions and applications.  SIAM review, 51(3):455–500, 2009.
  * [KBG18] Chiheon Kim, Afonso S Bandeira, and Michel X Goemans.  Stochastic block model for hypergraphs: Statistical limits and a semidefinite programming approach.  arXiv preprint arXiv:1807.02884, 2018.
  * [KDL80] Pieter M Kroonenberg and Jan De Leeuw.  Principal component analysis of three-mode data by means of alternating least squares algorithms.  Psychometrika, 45(1):69–97, 1980.
  * [Kol01] Tamara G Kolda.  Orthogonal tensor decompositions.  SIAM Journal on Matrix Analysis and Applications, 23(1):243–255, 2001.
  * [Kro08] Pieter M Kroonenberg.  Applied multiway data analysis, volume 702.  John Wiley & Sons, 2008.
  * [KS08] 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.
  * [KSS04] Amit Kumar, Yogish Sabharwal, and Sandeep Sen.  A simple linear time (1+/spl epsiv/)-approximation algorithm for k-means clustering in any dimensions.  In 45th Annual IEEE Symposium on Foundations of Computer Science, pages 454–462. IEEE, 2004.
  * [KSX19] Zheng Tracy Ke, Feng Shi, and Dong Xia.  Community detection for hypergraph networks via regularized tensor power iteration.  arXiv preprint arXiv:1909.06503, 2019.
  * [KX16] Vladimir Koltchinskii and Dong Xia.  Perturbation of linear forms of singular vectors under gaussian noise.  pages 397–423, 2016.
  * [LBW18] James Li, Jacob Bien, and Martin T Wells.  rtensor: An r package for multidimensional array (tensor) unfolding, multiplication, and decomposition.  Journal of Statistical Software, 87(10):1–31, 2018.
  * [LCL19] Jing Lei, Kehui Chen, and Brian Lynch.  Consistent community detection in multi-layer network data.  Biometrika, 2019.
  * [LL10] Nan Li and Baoxin Li.  Tensor completion for on-board compression of hyperspectral images.  In 2010 IEEE International Conference on Image Processing, pages 517–520. IEEE, 2010.
  * [LML+17] Thibault Lesieur, Léo Miolane, Marc Lelarge, Florent Krzakala, and Lenka Zdeborová.  Statistical and computational phase transitions in spiked tensor estimation.  In 2017 IEEE International Symposium on Information Theory (ISIT), pages 511–515. IEEE, 2017.
  * [LPV08] Haiping Lu, Konstantinos N Plataniotis, and Anastasios N Venetsanopoulos.  Mpca: Multilinear principal component analysis of tensor objects.  IEEE transactions on Neural Networks, 19(1):18–39, 2008.
  * [LR15] Jing Lei and Alessandro Rinaldo.  Consistency of spectral clustering in stochastic block models.  The Annals of Statistics, 43(1):215–237, 2015.
  * [LSF+14] Yuanyuan Liu, Fanhua Shang, Wei Fan, James Cheng, and Hong Cheng.  Generalized higher-order orthogonal iteration for tensor decomposition and completion.  In Advances in Neural Information Processing Systems, pages 1763–1771, 2014.
  * [LW20] Chanwoo Lee and Miaoyan Wang.  Tensor denoising and completion based on ordinal observations.  arXiv preprint arXiv:2002.06524, 2020.
  * [LXZL18] Xiaoshan Li, Da Xu, Hua Zhou, and Lexin Li.  Tucker tensor regression and neuroimaging analysis.  Statistics in Biosciences, pages 1–26, 2018.
  * [LYZ17] Tianqi Liu, Ming Yuan, and Hongyu Zhao.  Characterizing spatiotemporal transcriptome of human brain via low rank tensor decomposition.  arXiv preprint arXiv:1702.07449, 2017.
  * [LZ20a] Yuetian Luo and Anru R Zhang.  A schatten-q matrix perturbation theory via perturbation projection error bound.  arXiv preprint arXiv:2008.01312, 2020.
  * [LZ20b] Yuetian Luo and Anru R Zhang.  Tensor clustering with planted structures: Statistical optimality and computational limits.  arXiv preprint arXiv:2005.10743, 2020.
  * [MB18] Osman Asif Malik and Stephen Becker.  Low-rank tucker decomposition of large tensors using tensorsketch.  In Advances in Neural Information Processing Systems, pages 10096–10106, 2018.
  * [Mir60] Leon Mirsky.  Symmetric gauge functions and unitarily invariant norms.  The quarterly journal of mathematics, 11(1):50–59, 1960.
  * [MSC20] Xueyu Mao, Purnamrita Sarkar, and Deepayan Chakrabarti.  Estimating mixed memberships with sharp eigenvector deviations.  Journal of the American Statistical Association, (just-accepted):1–24, 2020.
  * [Ose11] Ivan V Oseledets.  Tensor-train decomposition.  SIAM Journal on Scientific Computing, 33(5):2295–2317, 2011.
  * [OT09] Ivan V Oseledets and Eugene E Tyrtyshnikov.  Breaking the curse of dimensionality, or how to use svd in many dimensions.  SIAM Journal on Scientific Computing, 31(5):3744–3759, 2009.
  * [OVW18] Sean O’Rourke, Van Vu, and Ke Wang.  Random perturbation of low rank matrices: Improving classical bounds.  Linear Algebra and its Applications, 540:26–59, 2018.
  * [PSB12] Evangelos E Papalexakis, Nicholas D Sidiropoulos, and Rasmus Bro.  From k-means to higher-way co-clustering: Multilinear decomposition with sparse latent factors.  IEEE transactions on signal processing, 61(2):493–506, 2012.
  * [PWB16] Amelia Perry, Alexander S Wein, and Afonso S Bandeira.  Statistical limits of spiked tensor models.  arXiv preprint arXiv:1612.07728, 2016.
  * [RCY11] Karl Rohe, Sourav Chatterjee, and Bin Yu.  Spectral clustering and the high-dimensional stochastic blockmodel.  The Annals of Statistics, 39(4):1878–1915, 2011.
  * [RM14] Emile Richard and Andrea Montanari.  A statistical model for tensor pca.  In Advances in Neural Information Processing Systems, pages 2897–2905, 2014.
  * [RYC19] Garvesh Raskutti, Ming Yuan, and Han Chen.  Convex regularization for high-dimensional multiresponse tensor regression.  The Annals of Statistics, 47(3):1554–1584, 2019.
  * [SDLF+17] Nicholas D Sidiropoulos, Lieven De Lathauwer, Xiao Fu, Kejun Huang, Evangelos E Papalexakis, and Christos Faloutsos.  Tensor decomposition for signal processing and machine learning.  IEEE Transactions on Signal Processing, 65(13):3551–3582, 2017\.
  * [SGL+19] Yiming Sun, Yang Guo, Charlene Luo, Joel Tropp, and Madeleine Udell.  Low-rank tucker approximation of a tensor from streaming data.  arXiv preprint arXiv:1904.10951, 2019.
  * [SL10] Berkant Savas and Lek-Heng Lim.  Quasi-newton methods on grassmannians and multilinear approximations of tensors.  SIAM Journal on Scientific Computing, 32(6):3352–3393, 2010.
  * [SN13] Andrey A Shabalin and Andrew B Nobel.  Reconstruction of a low-rank matrix in the presence of gaussian noise.  Journal of Multivariate Analysis, 118:67–76, 2013.
  * [SPL+09] Jimeng Sun, Spiros Papadimitriou, Ching-Yung Lin, Nan Cao, Shixia Liu, and Weihong Qian.  Multivis: Content-based social network exploration through multi-way visual analysis.  In Proceedings of the 2009 SIAM International Conference on Data Mining, pages 1064–1075. SIAM, 2009.
  * [SRT15] Parikshit Shah, Nikhil Rao, and Gongguo Tang.  Sparse and low-rank tensor decomposition.  In Advances in Neural Information Processing Systems, pages 2548–2556, 2015.
  * [SS07] Bernard N Sheehan and Yousef Saad.  Higher order orthogonal iteration of tensors (hooi) and its relation to pca and glram.  In Proceedings of the 2007 SIAM International Conference on Data Mining, pages 355–365. SIAM, 2007.
  * [Ste90] Gilbert W Stewart.  Matrix perturbation theory.  1990\.
  * [Ste98] Gilbert W Stewart.  Perturbation theory for the singular value decomposition.  Technical report, 1998.
  * [Ste06] Michael Stewart.  Perturbation of the svd in the presence of small singular values.  Linear algebra and its applications, 419(1):53–77, 2006.
  * [STFP12] Daniel L Sussman, Minh Tang, Donniell E Fishkind, and Carey E Priebe.  A consistent adjacency spectral embedding for stochastic blockmodel graphs.  Journal of the American Statistical Association, 107(499):1119–1128, 2012.
  * [SVBDL14] Laurent Sorber, Marc Van Barel, and Lieven De Lathauwer.  Tensorlab v2. 0.  Available online, January, 2014.
  * [Tso10] Charalampos E Tsourakakis.  Mach: Fast randomized tensor decompositions.  In Proceedings of the 2010 SIAM International Conference on Data Mining, pages 689–700. SIAM, 2010.
  * [Tuc66] Ledyard R Tucker.  Some mathematical notes on three-mode factor analysis.  Psychometrika, 31(3):279–311, 1966.
  * [Usc15] Andre Uschmajew.  A new convergence proof for the higher-order power method and generalizations.  Pacific Journal Of Optimization, 11(ARTICLE):309–321, 2015.
  * [Ver10] Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  arXiv preprint arXiv:1011.3027, 2010.
  * [VT03] M Alex O Vasilescu and Demetri Terzopoulos.  Multilinear subspace analysis of image ensembles.  In 2003 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2003. Proceedings., volume 2, pages II–93. IEEE, 2003.
  * [Vu11] Van Vu.  Singular vectors under random perturbation.  Random Structures & Algorithms, 39(4):526–538, 2011.
  * [VVM12] Nick Vannieuwenhoven, Raf Vandebril, and Karl Meerbergen.  A new truncation strategy for the higher-order singular value decomposition.  SIAM Journal on Scientific Computing, 34(2):A1027–A1052, 2012.
  * [Wan15] Rongrong Wang.  Singular vector perturbation under gaussian noise.  SIAM Journal on Matrix Analysis and Applications, 36(1):158–177, 2015.
  * [WBG16] Tao Wu, Austin R Benson, and David F Gleich.  General tensor spectral co-clustering for higher-order data.  In Advances in Neural Information Processing Systems, pages 2559–2567, 2016.
  * [WC14] Liqi Wang and Moody T Chu.  On the global convergence of the alternating least squares method for rank-one approximation to generic tensors.  SIAM Journal on Matrix Analysis and Applications, 35(3):1058–1072, 2014.
  * [WEAM19] Alexander S Wein, Ahmed El Alaoui, and Cristopher Moore.  The kikuchi hierarchy and tensor pca.  In 2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS), pages 1446–1468. IEEE, 2019.
  * [Wed72] PerAAke Wedin.  Perturbation bounds in connection with singular value decomposition.  BIT Numerical Mathematics, 12(1):99–111, 1972.
  * [Wey12] Hermann Weyl.  Das asymptotische verteilungsgesetz der eigenwerte linearer partieller differentialgleichungen (mit einer anwendung auf die theorie der hohlraumstrahlung).  Mathematische Annalen, 71(4):441–479, 1912.
  * [WLP15] 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.
  * [WZ19] Miaoyan Wang and Yuchen Zeng.  Multiway clustering via tensor block models.  In Advances in Neural Information Processing Systems, pages 713–723, 2019.
  * [Xu14] Yangyang Xu.  Fast algorithms for higher-order singular value decomposition from incomplete data.  arXiv preprint arXiv:1411.4324, 2014.
  * [Xu18] Yangyang Xu.  On the convergence of higher-order orthogonal iteration.  Linear and Multilinear Algebra, 66(11):2247–2265, 2018.
  * [XY17] Dong Xia and Ming Yuan.  On polynomial time methods for exact low-rank tensor completion.  Foundations of Computational Mathematics, pages 1–49, 2017.
  * [XYZ+05] Dong Xu, Shuicheng Yan, Lei Zhang, Hong-Jiang Zhang, Zhengkai Liu, and Heung-Yeung Shum.  Concurrent subspaces analysis.  In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), volume 2, pages 203–208. IEEE, 2005.
  * [XYZ17] Dong Xia, Ming Yuan, and Cun-Hui Zhang.  Statistically optimal and computationally efficient low rank tensor completion from noisy entries.  arXiv preprint arXiv:1711.04934, 2017.
  * [XZ19] Dong Xia and Fan Zhou.  The sup-norm perturbation of hosvd and low rank tensor denoising.  Journal of Machine Learning Research, 20(61):1–42, 2019.
  * [YWS14] Yi Yu, Tengyao Wang, and Richard J Samworth.  A useful variant of the davis–kahan theorem for statisticians.  Biometrika, 102(2):315–323, 2014.
  * [YXY+05] Shuicheng Yan, Dong Xu, Qiang Yang, Lei Zhang, Xiaoou Tang, and Hong-Jiang Zhang.  Discriminant analysis with tensor representation.  In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), volume 1, pages 526–532. IEEE, 2005.
  * [YXY+06] Shuicheng Yan, Dong Xu, Qiang Yang, Lei Zhang, Xiaoou Tang, and Hong-Jiang Zhang.  Multilinear discriminant analysis for face recognition.  IEEE Transactions on Image Processing, 16(1):212–220, 2006.
  * [YZ14] Ming Yuan and Cun-Hui Zhang.  On tensor completion via nuclear norm minimization.  Foundations of Computational Mathematics, pages 1–38, 2014.
  * [YZ16] Ming Yuan and Cun-Hui Zhang.  On tensor completion via nuclear norm minimization.  Foundations of Computational Mathematics, 16(4):1031–1068, 2016\.
  * [YZ17] Ming Yuan and Cun-Hui Zhang.  Incoherent tensor norms and their applications in higher order tensor completion.  IEEE Transactions on Information Theory, 63(10):6753–6766, 2017\.
  * [ZG01] Tong Zhang and Gene H Golub.  Rank-one approximation to high order tensors.  SIAM Journal on Matrix Analysis and Applications, 23(2):534–550, 2001.
  * [ZH18] Anru Zhang and Rungang Han.  Optimal sparse singular value decomposition for high-dimensional high-order data.  Journal of the American Statistical Association, page to appear, 2018.
  * [Zha19] Anru Zhang.  Cross: Efficient low-rank tensor completion.  The Annals of Statistics, 47(2):936–964, 2019.
  * [ZLRY20] Anru R Zhang, Yuetian Luo, Garvesh Raskutti, and Ming Yuan.  Islet: Fast and optimal low-rank tensor regression via importance sketching.  SIAM Journal on Mathematics of Data Science, 2(2):444–479, 2020\.
  * [ZLZ13] 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.
  * [ZT15] Qinqing Zheng and Ryota Tomioka.  Interpolating convex and non-convex tensor decompositions via the subspace norm.  In Advances in Neural Information Processing Systems, pages 3106–3113, 2015.
  * [ZWC+16] Hongya Zhao, Debby D Wang, Long Chen, Xinyu Liu, and Hong Yan.  Identifying multi-dimensional co-clusters in tensors based on hyperplane detection in singular vector spaces.  PloS one, 11(9), 2016.
  * [ZX18] Anru Zhang and Dong Xia.  Tensor svd: Statistical and computational limits.  IEEE Transactions on Information Theory, 64(11):7311–7338, 2018\.

Supplement to “A Sharp Blockwise Tensor Perturbation Bound for Orthogonal
Iteration”

Yuetian Luo, Garvesh Raskutti, Ming Yuan, and Anru R. Zhang

##  Appendix A Tensor Perturbation Bounds for HOOI in Asymmetric Case

In this section, we present a Corollary of Theorem 3 in the case where
$\Omega_{i}=\\{i\\},i=1,\ldots,d$, which appears frequently in practice.

######  Corollary 1 (Tensor Perturbation Bounds for HOOI in Asymmetric Case)

Consider the perturbation model (1) with
$\widetilde{{\mathbfcal{T}}},{\mathbfcal{T}},{\mathbfcal{Z}}\in\mathbb{R}^{p_{1}\times\cdots\times
p_{d}}$ and $\Omega_{i}=\\{i\\}$, $i=1,\ldots,d$. Define
$\mathcal{S}_{i}^{(-k)}:=\left\\{S\subseteq[d]\setminus\\{k\\}:|S|=i\right\\}$
as the set of all possible index sets with $i$ elements from
$[d]\setminus\\{k\\}$ and $\mathcal{S}_{0}^{(-k)}:=\emptyset$. For
$S\in\mathcal{S}_{i}^{(-k)}$, let $S^{c}=([d]\setminus\\{k\\})\setminus S$.
Now we define the blockwise errors as

$\begin{split}&\tau_{1}=\max_{k=1,\ldots,d}\tau_{1k},\quad\tau_{1k}=\left\|\left(\mathcal{M}_{k}({\mathbfcal{Z}}\times_{i\neq
k}{\mathbf{U}}_{i}^{\top})\right)_{\max(r_{k})}\right\|_{q},\quad
k=1,\ldots,d;\\\
&\tau_{j}=\max_{k=1,\ldots,d}\Big{\\{}\max_{S\in\mathcal{S}^{(-k)}_{j-1}}\sup_{\begin{subarray}{c}{\mathbf{V}}_{i}\in\mathbb{R}^{(p_{i}-r_{i})\times
r_{i}},\\\ \|{\mathbf{V}}_{i}\|_{q}\leq 1,i\in
S\end{subarray}}\left\|\left(\mathcal{M}_{k}({\mathbfcal{Z}}\times_{i\in
S}\left({\mathbf{U}}_{i\perp}{\mathbf{V}}_{i}\right)^{\top}\times_{i\in
S^{c}}{\mathbf{U}}_{i}^{\top}\right)_{\max(r_{k})}\right\|_{q}\Big{\\}},\\\
&\quad\text{ for }j=2,\ldots,d;\end{split}$

Denote the initialization error of
$\\{\widetilde{{\mathbf{U}}}_{k}^{(0)}\\}_{k=1}^{d}$ as
$e_{0}:=\max_{k=1,\ldots,d}\|\widetilde{{\mathbf{U}}}^{(0)\top}_{k\perp}{\mathbf{U}}_{k}\|_{q}$.
Assume the initialization error and the signal strength satisfy

$e_{0}\leq\frac{\sqrt{2}}{2}\,\text{ and }\,\lambda\geq
2^{\frac{d+4}{2}}\left(1+\frac{\sqrt{2}}{2}\right)^{d}\xi,$

where $\xi:=\underset{\|{\mathbfcal{Y}}\|_{{\rm HS}}\leq 1,{\rm
rank}({\mathbfcal{Y}})\leq(r_{1},\ldots,r_{d})}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\rangle$.

Then with inputs $\widetilde{{\mathbfcal{T}}}$,
$\\{\widetilde{{\mathbf{U}}}_{k}^{(0)}\\}_{k=1}^{d}$,
$\\{\Omega_{i}\\}_{i=1}^{d}$, the estimated mode-$k$ singular subspaces
updates in Algorithm 1 after $t$ iterations satisfy

$\max_{k\in[d]}\left\|\sin\Theta\left(\widetilde{{\mathbf{U}}}_{k}^{(t)},{\mathbf{U}}_{k}\right)\right\|_{q}\leq
2^{\frac{d+3}{2}}\frac{\tau_{1}}{\lambda}+\frac{e_{0}}{2^{t}}.$

Moreover, when $t_{\max}\geq\log(e_{0}\lambda/\tau_{1})\vee 1$, the outputs of
estimated mode-$k$ singular subspaces of Algorithm 1 satisfy

$\max_{k\in[d]}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}\leq\left(2^{\frac{d+3}{2}}+1\right)\frac{\tau_{1}}{\lambda},$
$\begin{split}\left\|\sin\Theta\left(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k}\right)\right\|_{q}\leq\frac{2}{\left(1-c^{*}(\tau_{1},\lambda,d)\right)^{\frac{d-1}{2}}}\left(\frac{\tau_{1k}}{\lambda}+\sum_{j=1}^{d-1}\frac{{d-1\choose
j}\left(2^{\frac{d+3}{2}}+1\right)^{j}\tau_{1}^{j}\tau_{j+1}}{\lambda^{j+1}}\right),\end{split}$

for $k=1,\ldots,d$, where
$c^{*}(\tau_{1},\lambda,d)=\left(2^{\frac{d+3}{2}}+2\right)^{2}\frac{\tau_{1}^{2}}{\lambda^{2}}\leq\frac{1}{2}$,
and the output of tensor reconstruction $\widehat{{\mathbfcal{T}}}$ satisfies

$\begin{split}\left\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\right\|_{{\rm
HS}}&\leq\left\|\llbracket{\mathbfcal{Z}};\widehat{{\mathbf{U}}}_{1}^{\top},\ldots,\widehat{{\mathbf{U}}}_{d}^{\top}\rrbracket\right\|_{{\rm
HS}}+\sum_{k=1}^{d}\left\|\widehat{{\mathbf{U}}}_{k\perp}^{\top}\mathcal{M}_{k}({\mathbfcal{T}})\right\|_{F}\\\
&\leq\left(1+2d\left(1-c^{*}(\tau_{1},\lambda,d)\right)^{-\frac{d-1}{2}}\right)\xi.\end{split}$

##  Appendix B Additional Proofs

###  B.1 Proof of Lemma 1

For generality, here we present the proof for order $d$ case. First we show
the first equivalent characterization.

$\begin{split}\underset{\|{\mathbfcal{Y}}\|_{{\rm HS}}\leq 1,{\rm
rank}({\mathbfcal{Y}})\leq(r_{1},\ldots,r_{d})}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\rangle&\geq\underset{\|{\mathbfcal{S}}\|_{{\rm
HS}}\leq
1,{\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}}}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{S}}\times_{1}{\mathbf{U}}_{1}\times\cdots\times_{d}{\mathbf{U}}_{d}\rangle\\\
&=\underset{\|{\mathbfcal{S}}\|_{{\rm HS}}\leq
1,{\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}}}{\sup}\langle{\mathbfcal{Z}}\times_{1}{\mathbf{U}}_{1}^{\top}\times\cdots\times{\mathbf{U}}_{d}^{\top},{\mathbfcal{S}}\rangle\\\
&=\underset{{\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}}}{\sup}\|{\mathbfcal{Z}}\times_{1}{\mathbf{U}}_{1}^{\top}\times\cdots\times_{d}{\mathbf{U}}_{d}^{\top}\|_{{\rm
HS}}.\end{split}$

On the other hand by Theorem 2 of [DLDMV00a], we have

${\mathbfcal{Y}}=\left({\mathbfcal{Y}}\times_{1}\widetilde{{\mathbf{U}}}_{1}^{\top}\times\cdots\times\widetilde{{\mathbf{U}}}_{d}^{\top}\right)\times_{1}\widetilde{{\mathbf{U}}}_{1}\times\cdots\times_{d}\widetilde{{\mathbf{U}}}_{d},$

where $\widetilde{{\mathbf{U}}}_{i}\in\mathbb{R}^{p_{i}\times r_{i}}$ is the
left singular space of $\mathcal{M}_{i}({\mathbfcal{Y}})$. Since
$\|{\mathbfcal{Y}}\times_{1}\widetilde{{\mathbf{U}}}_{1}^{\top}\times\cdots\times\widetilde{{\mathbf{U}}}_{d}^{\top}\|_{{\rm
HS}}\leq 1$, so

$\underset{\|{\mathbfcal{Y}}\|_{{\rm HS}}\leq 1,{\rm
rank}({\mathbfcal{Y}})\leq(r_{1},\ldots,r_{d})}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{Y}}\rangle\leq\underset{\|{\mathbfcal{S}}\|_{{\rm
HS}}\leq
1,{\mathbf{U}}_{i}\in\mathbb{O}_{p_{i},r_{i}}}{\sup}\langle{\mathbfcal{Z}},{\mathbfcal{S}}\times_{1}{\mathbf{U}}_{1}\times\cdots\times_{d}{\mathbf{U}}_{d}\rangle.$

So, we have proved the first equivalent characterization.

The second equivalent characterization comes directly from the consequence of
Theorem 4.1 and 4.2 of [DLDMV00b]. $\blacksquare$

###  B.2 Proof of Lemma 2

This proof idea of this lemma is to project ${\mathbfcal{Z}}$ onto orthogonal
subspaces ${\mathbf{U}}_{k}$ and ${\mathbf{U}}_{k\perp}$ at each modes.

$\begin{split}&{\mathbfcal{Z}}\times_{1}\widehat{{\mathbf{U}}}_{1}^{\top}\times\cdots\times\widehat{{\mathbf{U}}}_{d}^{\top}\\\
=&\left({\mathbfcal{Z}}\times_{1}(P_{{\mathbf{U}}_{1}}+P_{{\mathbf{U}}_{1\perp}})\times\cdots\times_{d}(P_{{\mathbf{U}}_{d}}+P_{{\mathbf{U}}_{d\perp}})\right)\times_{1}\widehat{{\mathbf{U}}}_{1}^{\top}\times\cdots\times_{d}\widehat{{\mathbf{U}}}_{d}^{\top}\\\
=&\left(\sum_{\Omega\subseteq[d]}{\mathbfcal{Z}}\times_{k\in\Omega}P_{{\mathbf{U}}_{k}}\times_{k\in\Omega^{c}}P_{{\mathbf{U}}_{k\perp}}\right)\times_{1}\widehat{{\mathbf{U}}}_{1}^{\top}\times\cdots\times_{d}\widehat{{\mathbf{U}}}_{d}^{\top}\\\
=&\sum_{\Omega\subseteq[d]}{\mathbfcal{Z}}\times_{k\in\Omega}\widetilde{{\mathbf{U}}}_{k}^{\top}P_{{\mathbf{U}}_{k}}\times_{k\in\Omega^{c}}\widetilde{{\mathbf{U}}}_{k}^{\top}P_{{\mathbf{U}}_{k\perp}}.\end{split}$

So by triangle inequality, we have

$\begin{split}\left\|\llbracket{\mathbfcal{Z}};\widehat{{\mathbf{U}}}_{1}^{\top},\ldots,\widehat{{\mathbf{U}}}_{d}^{\top}\rrbracket\right\|_{{\rm
HS}}&\leq\sum_{\Omega\subseteq[d]}\left\|{\mathbfcal{Z}}\times_{k\in\Omega}\widetilde{{\mathbf{U}}}_{k}^{\top}P_{{\mathbf{U}}_{k}}\times_{k\in\Omega^{c}}\widetilde{{\mathbf{U}}}_{k}^{\top}P_{{\mathbf{U}}_{k\perp}}\right\|_{{\rm
HS}}\\\
&\leq\sum_{\Omega\subseteq[d]}\left\|{\mathbfcal{Z}}\times_{k\in\Omega}{\mathbf{U}}_{k}^{\top}\times_{k\in\Omega^{c}}{\mathbf{U}}_{k\perp}^{\top}\right\|_{{\rm
HS}}\prod_{k\in\Omega}\|\widetilde{{\mathbf{U}}}_{k}^{\top}{\mathbf{U}}_{k}\|\prod_{k\in\Omega^{c}}\|\widetilde{{\mathbf{U}}}_{k}^{\top}{\mathbf{U}}_{k\perp}\|\\\
&\leq\sum_{\Omega\subseteq[d]}\left\|{\mathbfcal{Z}}\times_{k\in\Omega}{\mathbf{U}}_{k}^{\top}\times_{k\in\Omega^{c}}{\mathbf{U}}_{k\perp}^{\top}\right\|_{{\rm
HS}}\prod_{k\in\Omega^{c}}\|\sin\Theta(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k})\|\\\
&=\sum_{\Omega\subseteq\\{1,\ldots,d\\}}\theta_{\Omega}\prod_{k\in\Omega^{c}}\left\|\sin\Theta(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k})\right\|.\end{split}$

Here the second inequality is due to the fact that
$\|{\mathbfcal{Z}}\times_{i}{\mathbf{A}}{\mathbf{B}}\|_{{\rm
HS}}=\|{\mathbf{A}}{\mathbf{B}}\mathcal{M}_{i}({\mathbfcal{Z}})\|_{{\rm
HS}}\leq\|{\mathbf{A}}\|\|{\mathbf{B}}\mathcal{M}_{i}({\mathbfcal{Z}})\|_{{\rm
HS}}=\|{\mathbf{A}}\|\|{\mathbfcal{Z}}\times_{i}{\mathbf{B}}\|_{{\rm HS}}$ and
we apply iteratively for each mode with
${\mathbf{A}}=\widetilde{{\mathbf{U}}}_{k}^{\top}{\mathbf{U}}_{k}$ (or
$\widetilde{{\mathbf{U}}}_{k}^{\top}{\mathbf{U}}_{k\perp}$) and
${\mathbf{B}}={\mathbf{U}}_{k}^{\top}$ (or ${\mathbf{U}}_{k\perp}^{\top}$).
The third inequality is due to that
$\|\widetilde{{\mathbf{U}}}_{k}^{\top}{\mathbf{U}}_{k}\|\leq 1$ and
$\left\|\widetilde{{\mathbf{U}}}_{k}^{\top}{\mathbf{U}}_{k\perp}\right\|=\left\|\sin\Theta(\widehat{{\mathbf{U}}}_{k},{\mathbf{U}}_{k})\right\|$.
$\blacksquare$

###  B.3 Proof of Theorem 2

The proof is done by construction. Let’s denote
${\mathbfcal{I}}_{r}\in(\mathbb{R}^{r})^{\otimes d}$ as the order-$d$ identity
tensor with entries $(i,i,\ldots,i)$ to be $1$ and others are $0$. We
construct

${\mathbfcal{Z}}_{1}=\frac{\xi}{\sqrt{r}}{\mathbfcal{I}}_{r}\times_{1}\left(\begin{array}[]{c}{\mathbf{0}}_{r\times
r}\\\ {\mathbf{I}}_{r}\\\ {\mathbf{0}}_{(p_{1}-2r)\times
r}\end{array}\right)\times\cdots\times_{d}\left(\begin{array}[]{c}{\mathbf{0}}_{r\times
r}\\\ {\mathbf{I}}_{r}\\\ {\mathbf{0}}_{(p_{d}-2r)\times
r}\end{array}\right),$

where ${\mathbf{0}}_{m\times n}$ denotes a $m\times n$ matrix with all entries
to be $0$.

It is easy to check that $\underset{\|{\mathbfcal{Y}}\|_{{\rm HS}}\leq 1,{\rm
rank}({\mathbfcal{Y}})\leq(r,\ldots,r)}{\sup}\langle{\mathbfcal{Z}}_{1},{\mathbfcal{Y}}\rangle\leq\|{\mathbfcal{Z}}_{1}\|_{{\rm
HS}}=\xi$. Similarly we construct

${\mathbfcal{T}}_{1}=\frac{\xi}{\sqrt{r}}{\mathbfcal{I}}_{r}\times_{1}\left(\begin{array}[]{c}{\mathbf{I}}_{r}\\\
{\mathbf{0}}_{r\times r}\\\ {\mathbf{0}}_{(p_{1}-2r)\times
r}\end{array}\right)\times\cdots\times_{d}\left(\begin{array}[]{c}{\mathbf{I}}_{r}\\\
{\mathbf{0}}_{r\times r}\\\ {\mathbf{0}}_{(p_{d}-2r)\times
r}\end{array}\right).$

Also we let ${\mathbfcal{Z}}_{2}={\mathbfcal{T}}_{1}$ and
${\mathbfcal{T}}_{2}={\mathbfcal{Z}}_{1}$, and it is easy to check
$({\mathbfcal{T}}_{1},{\mathbfcal{Z}}_{1}),({\mathbfcal{T}}_{2},{\mathbfcal{Z}}_{2})\in\mathcal{F}_{r}(\xi)$.
At the same time, we have
${\mathbfcal{Z}}_{1}+{\mathbfcal{T}}_{1}={\mathbfcal{T}}_{2}+{\mathbfcal{Z}}_{2}$.
Thus

$\begin{split}\inf_{\widehat{{\mathbfcal{T}}}}\sup_{({\mathbfcal{T}},{\mathbfcal{Z}})\in\mathcal{F}_{r}(\xi)}\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm
HS}}&\geq\inf_{\widehat{{\mathbfcal{T}}}}\max\left\\{\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}_{1}\|_{{\rm
HS}},\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}_{2}\|_{{\rm HS}}\right\\}\\\
&\geq\frac{1}{2}\left(\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}_{1}\|_{{\rm
HS}}+\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}_{2}\|_{{\rm HS}}\right)\\\
&\geq\frac{1}{2}\|{\mathbfcal{T}}_{1}-{\mathbfcal{T}}_{2}\|_{{\rm
HS}}=\frac{\sqrt{2}}{2}\xi.\end{split}$

$\blacksquare$

###  B.4 Proof of Theorem 3

The proof is long and nontrivial. The main idea of the proof is to develop the
recursive error bound of $\widetilde{{\mathbf{U}}}_{k}^{(t+1)}$ i.e., the
estimate of ${\mathbf{U}}_{k}$ at iteration $t+1$, based on the error bound of
$\widetilde{{\mathbf{U}}}_{k}^{(t)}$, i.e., the estimate at iteration $t$. The
outline of the proof is the following: after a briefly introduction of
notations, the main proof could be divided into three steps.

  * •

Step 1: Recall in HOOI procedure, the update for the mode-$k$ singular
subspace satisfies

$\begin{split}\widetilde{{\mathbf{U}}}_{k}^{(t+1)}=&{\rm
SVD}_{r_{k}}\bigg{(}\mathcal{M}_{\bar{k}}\left({\mathbfcal{T}}\times_{i\in\underline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t+1)\top}_{i^{\prime}}\times_{\widecheck{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{k}\times_{i\in\overline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{i^{\prime}}\right)\\\
&\quad\quad\quad+\mathcal{M}_{\bar{k}}\left({\mathbfcal{Z}}\times_{i\in\underline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t+1)\top}_{i^{\prime}}\times_{\widecheck{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{k}\times_{i\in\overline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{i^{\prime}}\right)\bigg{)}.\end{split}$

To give an upper bound for $e_{t+1,k}$, we aim to give an upper bound for

$\left\|\left(\mathcal{M}_{\bar{k}}\left({\mathbfcal{Z}}\times_{i\in\underline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t+1)\top}_{i^{\prime}}\times_{\widecheck{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{k}\times_{i\in\overline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{i^{\prime}}\right)\right)_{\max(r_{k})}\right\|_{q}$

by using $\tau_{1},\ldots,\tau_{m}$, $e_{t},e_{t+1}$ in this step..

  * •

Step 2: After getting an upper bound for (20), in this step we use induction
to prove the following claim,

$e_{t}\leq 2^{(d+3)/2}\tau_{1}/\lambda+e_{0}/2^{t};t=0,1,\ldots.$

.

  * •

Step 3: Derive the error bound for
$\|\widehat{{\mathbfcal{T}}}-{\mathbfcal{T}}\|_{{\rm HS}}$ by the unified
quantity $\xi$.

For convenience, in this proof we denote

${\mathbf{T}}_{k}=\mathcal{M}_{\bar{k}}({\mathbfcal{T}}),\quad\widetilde{{\mathbf{T}}}_{k}=\mathcal{M}_{\bar{k}}(\widetilde{{\mathbfcal{T}}}),\quad{\mathbf{Z}}_{k}=\mathcal{M}_{\bar{k}}({\mathbfcal{Z}}),\quad
k=1,\ldots,m.$

Suppose

$e_{t}=\max_{k}e_{t,k},\quad
e_{t,k}=\left\|(\widetilde{{\mathbf{U}}}_{k\perp}^{(t)})^{\top}{\mathbf{U}}_{k}\right\|_{q},\quad
k=1,\ldots,m;t=0,1,\ldots.$ (31)

Step 1. Recall the procedure of HOOI that

$\begin{split}\widetilde{{\mathbf{U}}}_{k}^{(t+1)}=&{\rm
SVD}_{r_{k}}\left(\mathcal{M}_{\bar{k}}\left(\widetilde{{\mathbfcal{T}}}\times_{i\in\underline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t+1)\top}_{i^{\prime}}\times_{i\in\widecheck{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{k}\times_{i\in\overline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{i^{\prime}}\right)\right)\\\
=&{\rm
SVD}_{r_{k}}\bigg{(}\mathcal{M}_{\bar{k}}\left({\mathbfcal{T}}\times_{i\in\underline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t+1)\top}_{i^{\prime}}\times_{\widecheck{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{k}\times_{i\in\overline{\Omega}_{k}}\widetilde{{\mathbf{U}}}^{(t)\top}_{i^{\prime}}\right)\\\