Tensor SVD: Statistical and Computational Limits

∗

arXiv:1703.02724v3 [] 14 Apr 2018

Anru Zhang and Dong Xia
University of Wisconsin-Madison

Abstract
In this paper, we propose a general framework for tensor singular value decomposition
(tensor SVD), which focuses on the methodology and theory for extracting the hidden lowrank structure from high-dimensional tensor data. Comprehensive results are developed on
both the statistical and computational limits for tensor SVD. This problem exhibits three
different phases according to the signal-to-noise ratio (SNR). In particular, with strong SNR,
we show that the classical higher-order orthogonal iteration achieves the minimax optimal
rate of convergence in estimation; with weak SNR, the information-theoretical lower bound
implies that it is impossible to have consistent estimation in general; with moderate SNR,
we show that the non-convex maximum likelihood estimation provides optimal solution, but
with NP-hard computational cost; moreover, under the hardness hypothesis of hypergraphic
planted clique detection, there are no polynomial-time algorithms performing consistently in
general.

1

Introduction

There is no need to argue the importance of singular value decomposition (SVD) in data analysis.
As one of the most important tools in multivariate analysis, SVD along with the closely related
formulation, i.e., principal component analysis (PCA), have been a mainstay of data analysis
∗

Anru Zhang is Assistant Professor, Department of Statistics, University of Wisconsin-Madison, Madison,

WI 53706, E-mail: anruzhang@stat.wisc.edu; Dong Xia is Visiting Assistant Professor, Department of Statistics,
University of Wisconsin-Madison, Madison, WI 53706, E-mail: dongxia@stat.wisc.edu.

1

since more than a century ago, and widely used in various subjects. Attributed to the modern
high-dimensional data, the popularity of SVD and PCA continues to surge in the recent decades,
and many important variations, such as sparse SVD [1, 2, 3, 4], matrix denoising [5, 6, 7, 8], sparse
PCA [9, 10, 11], robust PCA [12], have been proposed and developed recently. Traditionally,
most of the SVD and PCA results focused on exploiting low-rank structures from datasets in
form of matrices.
Motivated by modern scientific research, tensors, or high-order arrays, have been actively
studied in machine learning, electrical engineering, and statistics. Some specific scientific applications involving tensor data include neuroimaging analysis [13, 14], recommender systems
[15, 16], computer vision [17, 18], topic modeling [19], community detection [20], hyperspectral
image compression [21], spatiotemporal gene expression [22], etc. A common objective in these
problems is to dig out the underlying high-order low-rank structure, such as the singular subspaces and the whole low-rank tensors, buried in the noisy observations. To achieve this goal, we
are in strong need of a statistical tool for tensor data that is the counterpart of regular singular
value decomposition for traditional order-2 datasets. Richard and Montanari [23], Hopkins et
al [24], Perry et al [25] considered a rank-1 spiked tensor SVD statistical model and proposed
various methods, including tensor unfolding and sum of square optimization (SOS). However,
as far as we know, the statistical framework for general rank-r high-order tensor SVD or PCA
was not well established or studied in the literature.
In this paper, we propose a general framework of tensor singular value decomposition (tensor
SVD). To be specific, suppose we are interested in a low-rank tensor X ∈ Rp1 ×p2 ×p3 , which is
observed with entry-wise corruptions as follows,
Y = X + Z.

(1)
iid

p1 ,p2 ,p3
Here Z is the p1 -by-p2 -by-p3 noisy tensor with {Zijk }i,j,k=1
∼ N (0, σ 2 ); X is a fixed tensor

with low Tucker ranks in the sense that all fibers of X along three directions (i.e., counterpart
of matrix columns and rows for tensors, also see Section 2.1 for formal definitions) lie in lowdimensional subspaces, say U1 , U2 , and U3 , respectively. Our goal is to estimate U1 , U2 , U3 , and
X from the noisy observation Y.
It is worth mentioning that the analog of this problem when X is an order-2 tensor, i.e.,
2

a matrix, has been previously studied in the context of matrix denoising in [5, 7, 8, 26]. For
the matrix denoising problem, the best low-rank matrix approximation provides the optimal
results, which can be calculated efficiently via singular value decomposition, as guaranteed by
the well-regarded Eckart-Young-Mirsky Theorem.
Although there have been significant efforts in developing methodologies and theories for
matrix SVD or matrix denoising, there is a paucity of literature on the analogous question for
tensors of order three or higher. In fact, SVD for high-order tensors is much more difficult than
its counterpart for matrices in various of aspects. First, tensors have more involved structures
along three or more ways, while the traditional tools for matrices could typically incorporate two
ways. As we will see later, one may achieve a sub-optimal result by simply ignoring the structure
beyond two ways. Second, many operations for matrices, such as operator norm, singular value
decomposition, are either not well defined or computational NP-hard for high-order tensors [27].
Third, high-order tensors often bring about high-dimensionality. For incidence, a 500-by-500-by500 tensor is comprised of more than 12,500,000 entries that impose significant computational
challenges. All these characteristics make the tensor SVD distinct from the classical matrix
problems.
The best low-rank tensor approximation, or equivalently the maximum likelihood estimation
(MLE), is a straightforward solution for tensor SVD. However, MLE is non-convex and computationally NP-hard in general (see, e.g. Hillar and Lim [27]). De Lathauwer, De Moor, and
Vandewalle instead introduced the higher order SVD (HOSVD) [28] and higher order orthogonal iteration (HOOI) [29], which aims at approximating the best low-rank approximation with
efficient spectral and power iteration method. Since then, HOSVD and HOOI have been widely
studied in the literature (see, e.g. [30, 31, 32, 33, 34]). However as far as we know, many basic
theoretical properties of these procedures, such as the error bound and the necessary iteration
times, still remain unclear.
In this paper, we develop comprehensive results on both the statistical and computational
limits for tensor SVD. To be specific, we establish upper bounds on estimation errors for
both higher-order orthogonal iteration (HOOI) and maximum likelihood estimator (MLE). It
is also shown that HOOI converges within a logarithm factor of iterations. Then the matching

3

information-theoretical lower bounds over a large class of low-rank tensors are correspondingly
introduced. To the best of our knowledge, we are among the first to develop the statistical
guarantees for both HOOI and MLE. Let the Tucker rank of X be (r1 , r2 , r3 ) (see formal definition in Section 3). The statistical and computational barriers of tensor SVD problem rely on
a key factor λ, i.e., the smallest non-zero singular values of matricizations of X (also see formal
definition in Section 3), which essentially measures the signal strength of the problem. When
p = min{p1 , p2 , p3 }, pk ≤ Cp, rk ≤ Cp1/2 for k = 1, 2, 3 and a constant C > 0, our main results
can be summarized into the following three phases according to signal-to-noise ratio (SNR): λ/σ.
1. When λ/σ = pα for α ≥ 3/4, the scenario is referred to as the strong SNR case. The
fast higher-order orthogonal iteration (HOOI) recovers U1 , U2 , U3 , and X with the minimax
optimal rate of convergence over a general class of low-rank tensors.
2. When λ/σ = pα for α < 1/2, we refer to this case as the weak SNR case, and propose
the minimax lower bound to show that there are no consistent estimators of U1 , U2 , U3 , or
X;
3. When λ/σ = pα for 1/2 ≤ α < 3/4, the scenario is referred to as the moderate SNR
case. We provide a computational lower bound to show that no polynomial time algorithm
can recover U1 , U2 , U3 consistently based on an assumption of hypergraphic planted clique
detection. Meanwhile, the maximum likelihood estimator, although being computational
intractable, achieves optimal rates of convergence over a general class of low-rank tensors.
It is also noteworthy that our results can be further generalized to fourth or higher order tensors,
or when the noise Z is i.i.d. sub-Gaussian distributed.
Our work is also related to several recent results in literature. For example, [23, 24, 35,
36] considered the extraction of rank-1 symmetric tensors from i.i.d. (symmetric) Gaussian
noise, which is a rank-1 special case of our tensor SVD model; [37, 38] considered the CP lowrank tensor decomposition based on noisy observations; [25] considered the statistical limit of
detecting and estimating a randomly sampled rank-one structure from a symmetric random
Gaussian tensor; [39, 40] considered the regularized tensor factorizations with/without sparsity;
[41] and [42] further considered non-negative tensor decomposition and robust tensor principal
4

component analysis; [22] focused on orthogonal decomposable tensor SVD problem; Lesieur et al
[43] considered a Bayesian symmetric spiked tensor estimation model – an approximate message
passing algorithm (AMP) was particularly introduced and the rigorous asymptotic analysis for
statistical and computational phase transitions were performed on high-order, symmetric, and
rank-1 tensor estimation. It should be noted that different from previous works, we perform
non-asymptotic analysis for tensor SVD, where the signal tensor X can be generally Tuckerrank-r, non-random, and asymmetric. Also, to the best of our knowledge, we are among the
first to provide a comprehensive analysis of both statistical and computational optimality of
tensor SVD.
The rest of the article is organized as follows. After a brief explanation for basic notations
and tensor algebra in Section 2.1, we state the fast higher-order orthogonal iteration and the
non-convex maximum likelihood estimation for tensor SVD in Section 2.2. The statistical limits
in the context of minimax optimality are provided for strong, weak, and moderate SNR cases
respectively in Section 3. Then we further discuss the computational barriers in moderate SNR
case in Section 4. Simulation results are provided in Section 5 to justify the theoretical results
of this paper. We briefly discuss the extension of the results to fourth or higher order tensors
and i.i.d. sub-Gaussian noise cases in Section 6. The proofs of all technical results are given in
Section 7 and the supplementary materials.

2
2.1

Tensor SVD: Methodology
Notations, Preliminaries, and Tensor Algebra

In this section, we start with basic notations and tensor algebra to be used throughout the
paper. For a, b ∈ R, let a ∧ b = min{a, b}, a ∨ b = max{a, b}. For two sequences {ai }, {bi },
if there are two constants C, c > 0 such that cai ≤ bi ≤ Cai for all i ≥ 1, we denote
a  b. C, c, C0 , c0 , . . . represent generic constants, whose actual values of these generic constants may vary from time to time. Particularly, the uppercase and lowercase letters represent
large and small constants, respectively. The matrices are denoted as capital letters, U1 , V1 , A,
etc. Especially, Op,r := {U ∈ Rp×r : U > U = Ir } is the set of all p-by-r matrices with or-

5

thonormal columns. For any matrix A ∈ Rp1 ×p2 , let σ1 (A) ≥ · · · ≥ σp1 ∧p2 (A) ≥ 0 be the
singular values in non-increasing order. We are particularly interested in the smallest singular value of A: σmin (A) = σp1 ∧p2 (A). In addition, the class of matrix Schatten q-norms
P
1/q
p1 ∧p2 q
will be used: kAkq =
σ
(A)
. Specific instances of Schatten q-norms include
j=1
j
qP
qP P
p1
p2
p1 ∧p2 2
2 =
the Frobenius norm (i.e., Schatten 2-norm), kAkF =
A
i=1
j=1 ij
j=1 σj (A),
and spectral norm (i.e., Schatten ∞-norm), kAk = σ1 (A) = maxv∈Rp2

kAvk2
kvk2 .

We also use

SVDr (A) to denote the leading r left singular vectors of A, so that SVDr (A) ∈ Op1 ,r . Define the projection operator PA = A(A> A)† A> . Here (·)† represents the psudo-inverse. If
A = U ΣV > is the SVD, PA can be equivalently written as PA = U U > . For any two matrices,
say U ∈ Rp1 ×r1 , V ∈ Rp2 ×r2 , we also let U ⊗ V ∈ R(p1 p2 )×(r1 r2 ) be their outer product matrix,
such that (U ⊗ V )[(i−1)p3 +j,(k−1)r3 +l] = Uik · Vjl , for i = 1, . . . , p1 , j = 1, . . . , p2 , k = 1, . . . , r1 , and
l = 1, . . . , r2 . We adopt R convention to represent submatrices: A[a:b,c:d] represents the a-to-b-th
rows, c-to-d-th columns of matrix A; we also use A[a:b,:] and A[:,c:d] to represent a-to-b-th full
rows of A and c-to-d-th full columns of A, respectively.
We use sin Θ distances to measure the difference between singular subspaces. To be specific,
for any two p × r matrices with orthonormal columns, say U and Û , we define the principal

angles between U and Û as Θ(U, Û ) = diag arccos(σ1 ), . . . , arccos(σr ) ∈ Rr×r , where σ1 ≥
. . . ≥ σr ≥ 0 are the singular values of U > Û . The Schatten q-sin Θ-norm is then defined as
!1/q
!1/q
r
r
X
X

q/2
k sin Θ(U, Û )kq =
=
sinq (arccos(σi ))
1 − σi2
, 1 ≤ q ≤ +∞.
i=1

i=1

The readers are referred to Lemma 3 in the supplementary materials and Lemma 1 in [26] for
more discussions on basic properties of sin Θ distances.
Throughout this paper, the boldface capital letters, e.g. X, Y, Z, note tensors. To simplify
the presentation, the main context of this paper is focused on third order tensor. The extension
to 4-th or higher tensors is briefly discussed in Section 6. The readers are also referred to [44]
for a more detailed tutorial of tensor algebra. For any tensor X ∈ Rp1 ×p2 ×p3 , define its mode-1
matricization as a p1 -by-(p2 p3 ) matrix M1 (X) such that
[M1 (X)]i,(j−1)p3 +k = Xijk ,

∀1 ≤ i ≤ p1 , 1 ≤ j ≤ p2 , 1 ≤ k ≤ p3 .

In other words, M1 (X) is composed of all mode-1 fibers, {(X[:,i2 ,i3 ] ) ∈ Rp1 : 1 ≤ i2 ≤ p2 , 1 ≤
i3 ≤ p3 }, of X. The mode-2 and mode-3 matricizations, i.e., M2 (X) ∈ Rp2 ×(p3 p1 ) and M3 ∈
6

Rp3 ×(p1 p2 ) , are defined in the same fashion. We also define the marginal multiplication ×1 :
Rp1 ×p2 ×p3 × Rr1 ×p1 → Rr1 ×p2 ×p3 as
X ×1 Y =

p1
X

!
Xi0 jk Yi,i0

i0 =1

.
1≤i≤r1 ,1≤j≤p2 ,1≤k≤p3

Marginal multiplications ×2 and ×3 can be defined similarly.
Different from matrices, there is no universal definition for tensor ranks. We particularly
introduce the following Tucker ranks (also called multilinear ranks) of X as
r1 = rank1 (X) = rank(M1 (X))
= dim(span{X[:,i2 ,i3 ] ∈ Rp1 : 1 ≤ i2 ≤ p2 , 1 ≤ i3 ≤ p3 }).
r2 = rank2 (X) and r3 = rank3 (X) can be similarly defined. Note that, in general, r1 , r2 , r3
satisfy r1 ≤ r2 r3 , r2 ≤ r3 r1 , r3 ≤ r1 r2 , but are not necessarily equal. We further denote rank(X)
as the triplet: (r1 , r2 , r3 ). The Tucker rank (r1 , r2 , r3 ) is also closely associated with the following
Tucker decomposition. Let U1 ∈ Op1 ,r1 , U2 ∈ Op2 ,r2 , U3 ∈ Op3 ,r3 be the left singular vectors of
M1 (X), M2 (X) and M3 (X) respectively, then there exists a core tensor S ∈ Rr1 ×r2 ×r3 such
that
X = S ×1 U1 ×2 U2 ×3 U3 ,

or

Xijk =

r1 X
r2 X
r3
X
i0 =1 j 0 =1

Si0 j 0 k0 (U1 )i,i0 (U2 )j,j 0 (U3 )k,k0 .

(2)

k0 =1

Expression (2) is widely referred to as the Tucker decomposition of X. Finally, to measure the
tensor estimation error, we introduce the following tensor Frobenius norm,
kXkF =

p1 X
p2 X
p3
X

2
Xijk

1/2

.

i=1 j=1 k=1

2.2

Maximum Likelihood Estimator and Higher-Order Orthogonal Iteration

In this section, we discuss the methodology for tensor SVD. Given the knowledge of Tucker
decomposition, the original tensor SVD model (1) can be cast as follows,
Y = X + Z = S ×1 U1 ×2 U2 ×3 U3 + Z,

i.i.d.

Z ∼ N (0, σ 2 ),

(3)

where U1 ∈ Op1 ,r1 , U2 ∈ Op2 ,r2 , U3 ∈ Op3 ,r3 , and S ∈ Rr1 ×r2 ×r3 . Our goal is to estimate
U1 , U2 , U3 , and X from Y. Clearly, the log-likelihood of Model (1) can be written (ignoring
7

the constants) as L (Y|X) = − σ12 kY − Xk2F , then it is straightforward to apply the maximum
likelihood estimator for estimation,
X̂mle =

arg min
rank(X)≤r1 ,r2 ,r3

Ûkmle = SVDrk (X̂mle ),

kY − XkF ,
k = 1, 2, 3.

Intuitively speaking, MLE seeks the best rank-(r1 , r2 , r3 ) approximation for Y in Frobenius
norm. By Theorems 4.1 and 4.2 in [29], MLE can be equivalently written as
Û1mle , Û2mle , Û3mle = arg max Y ×1 V1> ×2 V2> ×3 V3>
Vk ∈Opk ,rk

X̂

mle

2
F

,
(4)

= Y ×1 PÛ mle ×2 PÛ mle ×3 PÛ mle .
1

2

3

As we will illustrate later in Section 3, such estimators achieve optimal rate of convergence in
estimation errors. On the other hand, (4) is non-convex and computationally NP-hard even
when r = 1 (see, e.g., [27]). Then MLE may not be applicable in practice.
To overcome the computational difficulties of MLE, we consider a version of higher-order
orthogonal iteration (HOOI) [29]. The procedure includes three steps: spectral initialization,
power iteration and tensor projection. The first two steps produce optimal estimation of loadings
U1 , U2 , U3 . The final step yields an optimal estimator of the underlying low rank tensor X. It is
helpful for us to present the procedure of HOOI in details as follows.
Step 1 (Spectral initialization) Since U1 , U2 , and U3 respectively represent the singular subspaces
of M1 (X), M2 (X), and M3 (X), it is natural to perform singular value decomposition
(SVD) on Mk (Y) to obtain preliminary estimators for Uk :
(0)

Ûk
(0)

In fact, Ûk

= SVDrk (Mk (Y)) = the first rk left singular vectors of Mk (Y).

is exactly the higher-order SVD (HOSVD) estimator introduced by De Lath(0)

auwer, De Moor, and Vandewalle [28]. As we will show later, Ûk

serves as a good starting

point but not an optimal estimator for Uk .
Step 2 (Power Iteration) Then one applies power iterations to update the estimations. Given
(t−1)

Û2

(t−1)

, Û3

(t−1) >
) .

(Û3

(t−1) >
)

, Y can be denoised via mode-2 and 3 projections: Y ×2 (Û2

×3

As we will illustrate via theoretical analysis, the mode-1 singular subspace
8

of X is preserved while the amplitude of the noise is highly reduced after such the projection. Thus, for t = 1, 2, . . ., we calculate
(t)

(t−1) >

(t)

(t)

(t−1) >

(t)

(t)

(t)

Û1 = first r1 left singular vectors of M1 (Y ×2 (Û2

(t−1) >

) ×3 (Û3

Û2 = first r2 left singular vectors of M2 (Y ×1 (Û1 )> ×3 (Û3

) ),

) ),

(5)

Û3 = first r3 left singular vectors of M3 (Y ×1 (Û1 )> ×2 (Û2 )> ).
The iteration is stopped when either the increment is no more than the tolerance ε, i.e.,
(t)

(t)

(t)

Y ×1 (Û1 )> ×2 (Û2 )> ×3 (Û3 )>
− Y

(t−1) >
×1 (Û1
)

(t−1) >
×2 (Û2
)

F

(6)
(t−1) >

×3 (Û3

)

F

≤ ε,

or the maximum number of iterations is reached.
Step 3 (Projection) With the final estimates Û1 , Û2 , Û3 , it is natural to estimate S and X as
Ŝ = Y ×1 Û1> ×2 Û2> ×3 Û3> ,

X̂ = Ŝ ×1 Û1 ×2 Û2 ×3 Û3 = Y ×1 PÛ1 ×2 PÛ2 ×3 PÛ3 .

The procedure of HOOI is summarized in Algorithm 1. The further generalization to order-4 or
higher tensors SVD will be discussed in Section 6.

3

Statistical Limits: Minimax Upper and Lower Bounds

In this section, we develop the statistical limits for tensor SVD. Particularly, we analyze the
estimation error upper bounds of HOOI and MLE, then develop the corresponding lower bounds.
For any X ∈ Rp1 ×p2 ×p3 , denote λ = mink=1,2,3 σrk (Mk (X)) as the minimal singular values of
each matricization, which essentially measures the signal level in tensor SVD model. Suppose
the signal-to-noise ratio (SNR) is λ/σ = pα , where p = min{p1 , p2 , p3 }. Then the problem of
tensor SVD exhibits three distinct phases: α ≥ 3/4 (strong SNR), α < 1/2 (weak SNR), and
1/2 ≤ α < 3/4 (moderate SNR).
We first analyze the statistical performance of HOOI, i.e., Algorithm 1, under the strong
SNR setting that λ/σ ≥ Cp3/4 .
Theorem 1 (Upper Bound for HOOI). Suppose there exist constants C0 , c0 > 0 such that
pk ≤ C0 p, kXkF ≤ C0 σ exp(c0 p), rk ≤ C0 p1/2 , for p = min{p1 , p2 , p3 } and k = 1, 2, 3. Then
9

Algorithm 1 Higher Order Orthogonal Iteration (HOOI) [29]
1: Input: Y ∈ Rp1 ×p2 ×p3 , (r1 , r2 , r3 ), increment tolerance ε > 0, and maximum number of
iterations tmax .
2:

Let t = 0, initiate via matricization SVDs
(0)

Û1
3:
4:

= SVDr1 (M1 (Y)),

(0)

Û2

(0)

= SVDr2 (M2 (Y)),

Û3

= SVDr3 (M3 (Y)).

repeat
Let t = t + 1, calculate


(t)
(t−1) >
(t−1) >
Û1 = SVDr1 Mr1 (Y ×2 (Û2
) ×3 (Û3
) ) ,


(t)
(t)
(t−1) >
Û2 = SVDr2 Mr2 (Y ×1 (Û1 )> ×3 (Û3
) ) ,


(t)
(t)
(t)
Û3 = SVDr3 Mr3 (Y ×1 (Û1 )> ×2 (Û2 )> ) .

5:

until t = tmax or
(t)

(t)

(t)

Y ×1 (Û1 )> ×2 (Û2 )> ×3 (Û3 )>
− Y
6:

(t−1) >
×1 (Û1
)

(t−1) >
×2 (Û2
)

F
(t−1) >

×3 (Û3

)

F

Estimate and output:
(tmax )

Û1 = Û1

,

(tmax )

Û2 = Û2

,

X̂ = Y ×1 PÛ1 ×2 PÛ2 ×3 PÛ3 .

10

(tmax )

Û3 = Û3

;

≤ ε.

there exist absolute constants Cgap , C > 0, which does not depend on pk , rk , λ, σ, q, such that
whenever
λ/σ ≥ Cgap p3/4 ,

(i.e., in strong SNR case),

p
λ


∨ 1 iterations in Algorithm 1, the following upper bounds hold,
√


pk
−1/q
Erk
sin Θ Ûk , Uk
≤C
, k = 1, 2, 3, 1 ≤ q ≤ ∞,
(7)
λ/σ
q


2
kX̂ − Xk2F
(p1 + p2 + p3 ) ^
2
(8)
≤C
1 .
E X̂ − X ≤ Cσ (p1 r1 + p2 r2 + p3 r3 ) , E
λ2 /σ 2
F
kXk2F

after at most tmax = C log

Remark 1. In contrast to the error bound for final estimators Ûk in (7), an intermediate step
(0)

in the proof for Theorem 3 yields the following upper bound for initializations Ûk , i.e., the
output from Algorithm 1 Step 1,
√
−1/q
Erk

sin Θ



(0)
(0)
Ûk , Uk


q

≤C

pk
Cp3/2
+ 2 2,
λ/σ
λ /σ

k = 1, 2, 3.

(9)

Compared to Theorem 1, the bound in (9) is suboptimal as long as λ/σ = pα when 3/4 ≤ α < 1.
(0)

Thus, the higher-order SVD (HOSVD) Ûk

[28] may yield sub-optimal result. We will further

illustrate this phenomenon by numerical analysis in Section 5.
Remark 2. Especially when r = 1, Theorem 1 confirms the heuristic conjecture raised in
Richard and Montanari [23] that the tensor unfolding method yields reliable estimates for order3 spiked tensors provided that λ/σ > O(p3/4 ). Moreover, Theorem 1 further shows the power
iterations are necessary in order to refine the reliable estimates to minimax-optimal estimates.
Our result in Theorem 1 outperforms the ones by Sum-of-Squares (SOS) scheme (see, e.g.,
[24, 36]), where an additional logarithm factor on the assumption of λ is required. In addition,
the method we analyze here, i.e., HOOI, is efficient, easy to implement, and achieves optimal
rate of convergence for estimation error.
Remark 3. The strong SNR assumption (λ/σ ≥ Cp3/4 ) is crucial to guarantee the performance
of Algorithm 1. Actually, to ensure that Step 1 in Algorithm 1 provides meaningful initializations,
λ should be at least of order p3/4 according to our theoretical analysis.
Moreover, the estimators with high likelihood, such as MLE, achieve the following upper
bounds under weaker assumption that λ/σ ≥ Cp1/2 .
11

Theorem 2 (Upper Bound for Estimators with Large Likilihood and MLE). Suppose there exist
constants C0 , c0 > 0 such that pk ≤ C0 p, rk ≤ C0 p1/2 , kXkF ≤ C0 σ (exp(c0 p)), max{r1 , r2 , r3 } ≤
C0 min{r1 , r2 , r3 } for p = min{p1 , p2 , p3 } and k = 1, 2, 3. Suppose Ûk• ∈ Opk ,rk are estimators
satisfying
min kŶ − Ŝ• ×1 Û1• × Û2• ×3 Û3• k2F ≤ min kŶ − Ŝ ×1 U1 × U2 ×3 U3 k2F ,
Ŝ•

(10)

Ŝ

i.e., the likelihood value of Ûk• is no less than Uk . Then there exists a uniform constant Cgap > 0
(which does not depend on pk , rk , λ, σ, q) such that whenever
λ/σ ≥ Cgap p1/2 ,

(i.e., in moderate or strong SNR cases),

Û1• , Û2• , Û3• , and X̂• = Ŝ• ×1 Û1• ×2 Û2• ×3 Û3• satisfy
√
C pk
1/q
•
Erk
sin Θ(Ûk , Uk ) ≤
, k = 1, 2, 3,
λ/σ
q
E X̂• − X

2

F
kX̂• − Xk2F
E
kXk2F

1 ≤ q ≤ 2,

≤ Cσ 2 (p1 r1 + p2 r2 + p3 r3 ) ,


p1 + p2 + p3 ^
1 .
≤C
λ2 /σ 2

(11)

Especially, the upper bounds of (11) hold for maximum likelihood estimators (4).
Then we establish the lower bound for tensor SVD. We especially consider the following class
of general low-rank tensors,

Fp,r (λ) = X ∈ Rp1 ×p2 ×p3 : rankk (X) ≤ rk , σrk (Mk (X)) ≥ λ, k = 1, 2, 3 .

(12)

Here p = (p1 , p2 , p3 ), r = (r1 , r2 , r3 ) represent the dimension and rank triplets, λ is the smallest non-zero singular value for each matricization of X, which essentially measures the signal
strength of the problem. The following lower bound holds over Fp,r (λ).
Theorem 3 (Lower Bound). Suppose p = min{p1 , p2 , p3 }, max{p1 , p2 , p3 } ≤ C0 p, max{r1 , r2 , r3 } ≤
C0 min{r1 , r2 , r3 }, 4r1 ≤ r2 r3 , 4r2 ≤ r3 r1 , 4r3 ≤ r1 r2 , 1 ≤ rk ≤ pk /3, and λ > 0, then there exists
a universal constant c > 0 such that for 1 ≤ q ≤ ∞,
√ ^ 


pk
−1/q
1 ,
≥c
inf sup Erk
sin Θ Ũk , Uk
λ/σ
q
Ũk X∈Fp,r (λ)
inf

sup

X̂ X∈Fp,r (λ)

E X̂ − X

2
F

k = 1, 2, 3,

(13)

≥ cσ 2 (p1 r1 + p2 r2 + p3 r3 ) ,

kX̂ − Xk2F
≥c
inf sup E
kXk2F
X̂ X∈Fp,r (λ)
12




p1 + p2 + p3 ^
1 .
λ2 /σ 2

(14)

Remark 4. Theorem 3 contains two folds of meanings. First, when λ/σ ≤ cp1/2 for some
small constant c > 0, i.e., under weak SNR setting, the constant term dominates in (13), and
there are no consistent estimates for U1 , U2 , U3 . Secondly, when λ/σ ≥ Cp1/2 , i.e., under strong
√

and moderate SNR settings,

pk
λ/σ

dominates in (13), which provides non-trivial statistical lower

bounds for estimation errors.
We further define τ 2 = EkZk2F as the expected squared Frobenius norm of the whole tensor.
In summary, Theorems 1, 2, and 3 together yield the following statistical limits for tensor SVD.
1. Under strong SNR that λ/σ ≥ Cp3/4 (or λ/τ ≥ Cp−3/4 ), the higher-order orthogonal
iteration, i.e., Algorithm 1, provides minimax rate-optimal estimators for U1 , U2 , U3 , and
X.
√
pk
τ pk
 √
,

λ/σ
λ p1 p2 p3
√

inf

sup

Ûk X∈Fp,r (λ)

−1/q
Erk

sin Θ(Ûk , Uk )

q

k = 1, 2, 3, 1 ≤ q ≤ +∞,

τ 2 (p1 r1 + p2 r2 + p3 r3 )
,
p1 p2 p3
F
X̂ X∈Fp,r (λ)
  2


kX̂ − Xk2F
τ (p1 + p2 + p3 )
p1 + p2 + p3
inf sup E
∧1 

∧1 .
λ2 /σ 2
λ2 p1 p2 p3
kXk2F
X̂ X∈Fp,r (λ)
inf

sup

E X̂ − X

2

 σ 2 (p1 r1 + p2 r2 + p3 r3 ) 

(15)
2. Under moderate SNR that Cp1/2 ≤ λ/σ ≤ cp3/4 (or Cp−1 ≤ λ/τ ≤ cp−3/4 ), the estimators
with high likelihood (10), including the MLE (4), are minimax rate-optimal. The rate here
is exactly the same as (15).
3. Under weak SNR that λ/σ ≤ cp1/2 (or λ/τ ≤ cp−1 ), there are no consistent estimators for
U1 , U2 , U3 , or X.
However, as we have discussed in Section 2.2, MLE is not applicable even with moderate
dimension. It is still crucial to know whether there is any fast and efficient algorithm for tensor
SVD under moderate SNR setting.

13

4

Computational Limits in Moderate SNR Case

In this section, we focus on the computational aspect of tensor SVD under moderate SNR setting.
If λ/σ = pα with p = min{p1 , p2 , p3 }, α < 3/4, we particularly develop the computational lower
bound to show that every polynomial-time algorithm is statistically inconsistent in estimating U1 ,
U2 , U3 , and X based on computational hardness assumption. In recent literature, we have seen
achievements in obtaining computational lower bounds via computational hardness assumptions
for many problems, such as sparse PCA [45, 46, 47, 48], submatrix localization [49, 50, 51],
tensor completion [52], sparse CCA [48], and community detection [53]. The computational
hardness assumptions, such as planted clique detection and Boolean satisfiability, has been
widely studied and conjectured that no polynomial-time algorithm exists under certain settings.
For tensor SVD, our computational lower bound is established upon the hardness hypothesis of
hypergraphic planted clique detection, which is discussed in details in the next section.

4.1

Planted clique detection in hypergraphs

Let G = (V, E) be a graph, where V = {1, 2, . . . , N } and E are the vertex and edge sets,
respectively. For a standard graph, the edge e = (i, j) ∈ E indicates certain relation exists
between vertices i and j in V . A 3-hypergraph (or simply noted as hypergraph, without causing
confusion) is a natural extension, where each hyper-edge is represented by an unordered group
of three different vertices, say e = (i, j, k) ∈ E. Given a hypergraph G = (V, E) with |V | = N ,
its adjacency tensor A ∈ {0, 1}N ×N ×N is defined as

Aijk =



1, if (i, j, k) ∈ E;

0, otherwise.

We denote the Erdős-Rényi hypergraph of N vertices as G3 (N, 1/2), if for each 1 ≤ i < j < k ≤
N , (i, j, k) is included into the hyper-edge set independently with probability 1/2. For V1 ⊂ V
and certain integer 1 ≤ κN ≤ |V1 |, we use G3 (N, 1/2, κN , V1 ) to denote a random hypergraph
where a clique of size κN is planted inside V1 . More precisely, we first sample a random graph
from G3 (N, 1/2), then pick κN vertices uniformly at random from V1 , denote them as C, and
connecting all hyper-edges (i, j, k) for all distinct triplets i, j, k ∈ C. Conventionally, the planted
14

clique detection is referred to as the problem for distinguishing whether there is any planted
clique hidden in the Erdős-Rényi graph. To simplify our analysis in tensor SVD later, we propose
a slightly different version of hypergraphic planted clique detection problem as follows.
Definition 1. Let G be drawn from either G3 (N, 1/2, κN , V1 ) or G3 (N, 1/2, κN , V2 ), where V1 =

{1, 2, . . . , bN/2c} and V2 = bN/2c + 1, bN/2c + 2, . . . , N . The hypergraphic planted clique
detection problem, noted as PC3 (N, κN ), refers to the hypothesis testing problem
H0 : G ∼ G3 (N, 1/2, κN , V1 )

vs.

H1 : G ∼ G3 (N, 1/2, κN , V2 ).

(16)

Given a hypergraph G sampled from either H0 or H1 with adjacency tensor A ∈ {0, 1}N ×N ×N ,
let ψ(·) : {0, 1}N ×N ×N 7→ {0, 1} be a binary-valued function on A such that ψ(A) = 1 indicates
rejection of H0 . Then the risk of test ψ is defined as the sum of Type-I and II errors,


RN,κN (ψ) = PH0 ψ(A) = 1 + PH1 ψ(A) = 0 .
Put it differently, given a random hypergraph G ∼ H0 or H1 , our goal is to identify whether
the clique is planted in the first or second half of vertices.
When we replace the hyper-edges (involving three vertices each) of G3 (N, 1/2, κN , V1 ) by
the regular edges (involving two vertices each), the above hypergraphic planted clique detection
becomes the traditional planted clique detection problem. To provide circumstantial evidence
to the hardness of PC3 (N, κN ), it is helpful for us to review some well-known results of the traditional planted clique detection here. First, the difficulty of traditional planted clique detection

depends crucially on the planted clique size: κN . [54] and [55] showed that if κN = o log N ,
it is statistically impossible to determine whether a planted clique exists since a random graph
√
G ∼ G2 (N, 1/2) contains a clique of size 2 log N with high probability. When κN ≥ C N , it has
been shown that the planted clique can be located by performing polynomial-time operations
√
by spectral methods [56, 57]. If the size clique further increases, say κN ≥ C N log N , [58]
developed an algorithm to find exactly the planted clique with high probability in polynomial
√
time. However, when log N  κN  N , there is still no known polynomial-time algorithm
for planted clique detection, and it is currently widely conjectured by the theoretical computer
science and graph theory community that such polynomial-time algorithm may not exist (see
[59], [60],[55] and references therein).
15

When moving to hypergraphs, the hardness of PC3 (N, κN ), to the best of our knowledge,
remains unclear. In an extreme case of exhaustive search, it needs an exponential number of op
erations, i.e., κNN , to verify a solution. In addition, the performance of the simple matricizationspectral method (which shares similar idea as the proposed Algorithm 1) highly depends on the
size of the clique κN . We particularly have the following Proposition 1.
Proposition 1. Suppose G ∼ G3 (N, 1/2, κN , V1 ), so there exists C ⊆ V1 as a planted clique
of size κN with uniform random position. Let A be the corresponding adjacency tensor, and
1C ∈ R|V1 | be the indicator for the hidden clique that (1C )i = 1{i∈C} . We further partition
V1 = {1, . . . , bN/2c} into three equal subsets: Dk = {bkN/6c+1, . . . , b(k+1)N/6c} for k = 1, 2, 3.
Then we can calculate ûk ∈ R|V1 | as the leading left singular vector of Mk (2 · A[D1 ,D2 ,D3 ] −
1|D1 |×|D2 |×|D3 | ), where 1|D1 |×|D2 |×|D3 | is a |D1 |-by-|D2 |-by-|D3 | tensor with all entries 1. If the
sequence {κN } satisfies lim inf N →∞

κN
N 1/2

= ∞, then
d

sin Θ (ûk , (1C )Dk ) → 0,

as N → ∞,

k = 1, 2, 3.

In another word, the angle between ûk and (1C )Dk tends to 0 in probability.
Remark 5. For technical convenience, we partition V1 into three parts and perform SVD
on Mk (2A[D1 ,D2 ,D3 ] − 1|D1 |×|D2 |×|D3 | ) to ensure that most of the entries of Mk (A) are i.i.d.
Rademacher distributed.
Proposition 1 suggests that ûk can be used to locate C when κN  N 1/2 . However, the
theoretical analysis in Proposition 1 fails when κN = N (1−τ )/2 for τ > 0, and we conjecture that
such computational barrier is essential. Particularly, we propose the following computational
hardness assumption on hypergraphic planted clique detection.
log√
κN
≤ (1 − τ ) and any
Hypothesis 1. H(τ ). For any sequence {κN } such that lim sup log
N
N →∞

sequence of polynomial-time tests {ψN },
1
lim inf RN,κN (ψN ) ≥ .
N →∞
2

16

4.2

The computational lower bound of tensor SVD

Now we are ready to develop the computational lower bound for tensor SVD based on Hypothesis
H(τ ). Recall
Y = X + Z ∈ Rp1 ×p2 ×p3 ,

X = S ×1 U1 ×2 U2 ×3 U3 ,

iid

Z ∼ N (0, σ 2 ).

To better present the asymptotic argument, we add a superscript of dimension, p = min{p1 , p2 , p3 },
(p)

to the estimators, i.e., Ûk , X̂(p) . The computational lower bound is then presented as below.
Theorem 4 (Computational Lower Bound). Suppose hypergraphic planted clique assumption
H(τ ) holds for some τ ∈ (0, 1). Then there exist absolute constants c0 , c1 > 0 such that if
 3(1−τ )/4 
(p)
λ/σ ≤ c0 p√log 3p , for any integers r1 , r2 , r3 ≥ 1 and any polynomial time estimators Ûk ,
X̂(p) , the following inequalities hold
lim inf

sup

p→∞ X∈F (λ)
p,r

lim inf

(p)

E sin Θ Ûk , Uk

sup

p→∞ X∈F (λ)
p,r



2

≥ c1 ,

k = 1, 2, 3,

(17)

EkX̂(p) − Xk2F
≥ c1 .
kXk2F

(18)

Remark 6. For technical reasons, there is an additional logarithmic factor in the condition
 3(1−τ )/4 
λ/σ ≤ c0 p√log 3p , compared with the statistical lower bound in Theorem 3. Since τ is a strictly
positive number, the effect of logarithmic factor is dominated by pc for any c > 0 asymptotically.
Theorem 4 illustrates the computational hardness for tensor SVD under moderate scenario
that λ/σ = pα , 1/2 ≤ α < 3/4, if the hypergraphic planted clique assumption H(τ ) holds for
any τ > 0.

5

Simulations

In this section, we further illustrate the statistical and computational limits for tensor SVD via
numerical studies.
(0)

We first consider the average Shatten q-sin Θ-norm losses for initial estimators Ûk (HOSVD)
and final estimators Ûk (HOOI) under the following simulation setting. For given (p, r, λ), we
17

(p, r, λ)

l1 (Û )

l1 (Û (0) )

l2 (Û )

l2 (Û (0) )

l5 (Û )

l5 (Û (0) )

l∞ (Û )

l∞ (Û (0) )

(50, 5, 20)

1.1094

2.1192

0.5194

1.0535

0.3572

0.7991

0.3286

0.7699

(50, 5, 50)

0.4297

0.5243

0.2016

0.2519

0.1392

0.1815

0.1283

0.1713

(50, 10, 20)

2.4529

4.5208

0.8179

1.5674

0.4629

0.9611

0.3955

0.8762

(50, 10, 50)

0.9111

1.1210

0.3030

0.3771

0.1707

0.2175

0.1452

0.1890

(100, 5, 40)

0.7952

1.5649

0.3695

0.7707

0.2509

0.5778

0.2294

0.5543

(100, 5, 60)

0.5301

0.8132

0.2463

0.3938

0.1673

0.2878

0.1530

0.2731

(100, 10, 40)

1.7448

3.5371

0.5688

1.1943

0.3087

0.7015

0.2554

0.6246

(100, 10, 60)

1.1466

1.8055

0.3735

0.6015

0.2021

0.3427

0.1660

0.2975

Table 1: The average Schatten q-sin Θ loss of the final estimations Ûk and the spectral ini(0)

tializations Ûk based on 100 repetitions. Here, p1 = p2 = p3 = p, r1 = r2 = r3 = r,
P
lq (Û ) = 13 3k=1 k sin Θ(Ûk , Uk )kq .
let p = p1 = p2 = p3 , r = r1 = r2 = r3 , generate Ũk ∈ Rpk ×rk as i.i.d. standard Gaussian
matrices, and apply QR decomposition on Ũk and assign the Q part to Uk . In other words,
the singular subspaces U1 , U2 , U3 are drawn randomly from Haar measure. Then we construct
S̃ ∈ Rr1 ×r2 ×r3 as an i.i.d. Gaussian tensor, and rescale as S = S̃ ·

λ
mink=1,2,3 σrk (Mk (S̃))

to ensure

mink=1,2,3 σrk (Mk (X)) ≥ λ. Next, we construct Y = X + Z, where the signal tensor X =
S ×1 U1 ×2 U2 ×3 U3 , the noise tensor Z are drawn from i.i.d. standard Gaussian distribution.
We apply Algorithm 1 to Y, and record the average numerical performance under various values
of p, r, λ. The results based on 100 replications are shown in Table 1. We can clearly see that
the power iterations (Step 2 in Algorithm 1, i.e., HOOI) significantly improve upon spectral
initializations (Step 1 in Algorithm 1, i.e., HOSVD) in different Shatten-q sin Θ losses under
various settings.
Then we consider another setting that X has three different dimensions. Specifically, we
generate Y = X + Z by the same scheme as the previous setting with varying (p1 , p2 , p3 ) and
fixed r1 = r2 = r3 = 5. We repeat the experiment for 100 times, then record the average
estimation errors in Table 2. Again, we can see HOOI performs well under various values of
dimensions.

18

(p1 , p2 , p3 , λ)

l∞ (Û1 )

l2 (Û1 )

l∞ (Û2 )

l2 (Û2 )

l∞ (Û3 )

l2 (Û3 )

kX̂ − XkF

kX̂−XkF
kXkF

(20, 30, 50, 20)

0.2082

0.3032

0.2530

0.3858

0.3109

0.4975

24.7037

0.3276

(20, 30, 50, 100)

0.0409

0.0596

0.0498

0.0761

0.0641

0.1017

23.5708

0.0631

(30, 50, 100, 20)

0.2674

0.4036

0.3354

0.5247

0.4456

0.7252

33.6219

0.4479

(30, 50, 100, 100)

0.0490

0.0753

0.0640

0.1012

0.0911

0.1469

30.9540

0.0822

(100, 200, 300, 50)

0.1840

0.2982

0.2551

0.4301

0.3161

0.5155

57.8482

0.3090

(100, 200, 300, 100)

0.0940

0.1506

0.1259

0.2117

0.1638

0.2627

55.9009

0.1505

(200, 300, 400, 50)

0.2579

0.4335

0.3331

0.5523

0.3420

0.6017

72.2912

0.4026

(200, 300, 400, 150)

0.0825

0.1389

0.1076

0.1739

0.1277

0.2024

68.0305

0.1199

Table 2: The average spectral and Frobenius sin Θ loss for Û1 , Û2 , Û3 and average Frobenius
loss for X̂ under various settings. Here l∞ (Ûk ) = k sin Θ(Ûk , Uk )k, l2 (Û ) = k sin Θ(Ûk , Uk )kF .
Next, we illustrate the phase transition phenomenon of tensor SVD. Let X = S×1 U1 ×2 U2 ×3
U3 be a p-by-p-by-p tensor, where U1 , U2 , U3 are randomly generated p-by-r orthogonal matrices
from Haar measure and S ∈ Rr×r×r is a fixed diagonal tensor such that Si,j,k = pα · 1{i=j=k} , 1 ≤
i, j, k ≤ r for α ∈ [0.4, 0.9]. Then pα is the signal strength in our context. The entries of Z
√ √
are generated as either i.i.d. N (0, 1) or Unif[− 3, 3], which are sub-Gaussian, mean 0, and
variance 1. To demonstrate the phase transitions at both p3/4 and p1/2 , ideally one wishes to
implement both MLE and HOOI. Since MLE, i.e., the best low-rank approximation estimator
(4), is computationally intractable, we instead consider the following oracle warm-start HOOI
to obtain an approximation for MLE: suppose an oracle provides a warm start as
(0)warm

Ûk

1
1
= √ Uk + √ Uk0 ,
2
2

k = 1, 2, 3,

where Uk is the true underlying loading and Uk0 is a p-by-r random orthonormal matrix in the
complementary space of Uk . {Uk0 }3k=1 here are generated based on the following scheme: first
calculate Uk⊥ ∈ Op,p−r as the orthogonal complement of Uk , then construct Uk0 = Uk⊥ O for
some random orthogonal matrix O ∈ Op−r,r . Based on the oracle warm-start, we apply Steps 2
and 3 of Algorithm 1 to obtain the warm-start HOOI estimator Ûkwarm as an approximation for
MLE.
We let p vary from 50 to 100, r = 5, and apply both the spectral-start HOOI (i.e., the
original HOOI and Algorithm 1) and the oracle warm-start HOOI. The average spectral sin Θ
19

loss, i.e., l∞ (Û ) =

1
3

P3

k=1 k sin Θ(Ûk , Uk )k

from 100 repetitions are presented in Figure 1, where

the upper panel and lower panel correspond to the i.i.d. Gaussian noise and i.i.d. uniform noise
cases, respectively. Both panels of Figure 1 clearly demonstrate the phase transition effects: the
estimation error significantly decreases around SNR = p3/4 and around SNR = p1/2 for spectralstart HOOI and oracle warm-start HOOI, respectively. This exactly matches our theoretical
findings in Section 3. In addition, there is little difference between two plots in the upper and
lower panels, which implies that the statistical estimation error for tensor SVD mainly relies on
the SNR and is less influenced by the particular sub-Gaussian noise type.

6

Discussions: Further Generalizations

In this article, we propose a general framework for tensor singular value decomposition (tensor
SVD), which focuses on extracting the underlying Tucker low-rank structure from the noisy
tensor observations. We provide a comprehensive analysis on tensor SVD in aspects of both
statistics and computation. The problem exhibits three distinct phases according to the signalto-noise ratio (SNR): with strong SNR, the higher order orthogonal iteration (HOOI) performs
efficiently and achieves statistical optimal results; with weak SNR, no method performs consistently; with moderate SNR, the estimators with high likelihood, such as the computational
intractable MLE, perform optimally in statistical convergence rate, and no polynomial algorithm can do so unless we have a polynomial-time algorithm for the hypergraphic planted clique
problem.
The results of this paper are mainly presented under i.i.d. Gaussian noise. However, when
the noise is more generally i.i.d. sub-Gaussian distributed, say
iid

Z ∼ Z,

where

kZkψ2 = sup q −1/2 (E|Z|q )1/q ≤ σ,
q≥1

we can derive the upper bound results similarly as Theorems 1 and 2, as the proofs of main
technical tools, including Lemmas 5 and 8, still hold for i.i.d. sub-Gaussian noise case.
We have also focused our presentations mainly on order-3 tensors throughout this article.
The results can be additionally generalized to order-d tensor SVD for any d ≥ 2. Suppose one

20

0.9
0.8

l∞ (Û )

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.4

warm start: p=50
warm start: p=80
warm start: p=100
spectral start: p=50
spectral start: p=80
spectral start: p=100
0.5

0.6

0.7

0.8

0.9

0.7

0.8

0.9

α

0.9
0.8

l∞ (Û )

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.4

warm start: p=50
warm start: p=80
warm start: p=100
spectral start: p=50
spectral start: p=80
spectral start: p=100
0.5

0.6

α
Figure 1: Phase transitions in tensor SVD at SNR = p.5 and p.75 . Upper panel: Gaussian noise
√ √
N (0, 1). Lower panel: uniform noise Unif[− 3, 3].

21

observes an order-d tensor as follows,
Y = X + Z,

X = S ×1 U1 · · · ×d Ud ,

(19)

where Y, X, Z ∈ Rp1 ×···×pd , Uk ∈ Opk ,rk for k = 1, . . . , d, and S ∈ Rr1 ×···×rd . The higher-order
orthogonal iteration can be similarly written as follows [29],
Step 1 Initialize by singular value decomposition of each matricizations,
(0)

Ûk

= SVDrk (Mk (Y)) ,

k = 1, . . . , d.

Step 2 Let t = 0, 1, . . ., update the estimates for Uk sequentially for k = 1, . . . , d,

(t+1)
(t+1) >
(t+1)
Ûk
= SVDrk Y ×1 (Û1
) × · · · ×(k−1) (Ûk−1 )>

(t) >
(t) >
×(k+1) (Ûk+1 ) × · · · ×d (Ûd ) .
The iteration is continued until convergence or maximum number of iteration is reached.
Step 3 With the final estimators {Ûk }dk=1 from Step 2, one estimates X as
X̂ = Y ×1 PÛ1 × · · · ×d PÛd .
Meanwhile, the non-convex maximum likelihood estimates can be written as
(Û1mle , . . . , Ûdmle ) = arg max Y ×1 V1> × · · · ×d Vd>
Vk ∈Opk ,rk
k=1,...,d

F

,
(20)

X̂mle =Y ×1 PÛ mle × · · · ×d PÛ mle .
1

d

Again, let λ = min1≤k≤d σrk (Mk (X)) measure the signal strength. For fixed d, when p =
min{p1 , . . . , pd }, max{p1 , . . . , pd } ≤ Cp, rk ≤ Cp1/(d−1) , similarly as the proofs for Theorems 1,
2, and 3, it is possible to show under strong SNR case, where λ/σ = pα for α ≥ d/4, HOOI
achieves optimal rate of convergence over the following class of low-rank tensors

Fp,r (λ) = X ∈ Rp1 ×···×pd : σrk (Mk (X)) ≥ λ, k = 1, . . . , d ;
under the weak SNR where λ/σ = pα for α < 1/2, it is impossible to generally have consistent
estimators for U1 , . . . , Ud , or X; under the moderate SNR, where λ/σ = pα for
22

1
2

≤ α <

d
4,

the estimators with high likelihood, such as MLE, achieves optimal statistical performance;
while one can develop computational lower bound with computational hardness assumption of a
higher-order hypergraphic planted clique detection problem similarly as Theorem 4. We can also
see that the gap between statistical and computational limits vanishes if d = 2. This coincides
with the previous results in matrix denoising literature (see, e.g. [7, 8, 26]), where standard
singular value decomposition achieves both statistical optimality and computational efficiency.
Additionally, if d grows rather than stays as a fixed constant, the asymptotics of tensor SVD
in both statistical and computational aspects will be an interesting future project.

7

Proofs

We collect the proofs in this section for the main results in this paper. To be specific, the proof
for Theorems 1, 3, and 4 will be presented in Sections 7.1, 7.3, and 7.2, respectively. Proofs for
Theorem 2, Proposition 1, and additional technical lemmas are postponed to the supplementary
materials.

7.1

Proof of Theorem 1

We first consider the proof for Theorem 1. Throughout the proof, we assume the noise level
σ 2 = 1 without loss of generality. For convenience, we denote
X1 = M1 (X),

X2 = M2 (X),

X3 = M3 (X)

as the matricizations of X. We also denote Y1 , Y2 , Y3 , Z1 , Z2 , Z3 in the similar fashion. We also
let r = max{r1 , r2 , r3 }.
We divide the proof into steps.
1. In this first step, we consider the performance of initialization step, we particularly prove
that for any small constant c0 > 0, there exists large constant Cgap > 0 such that whenever
λ ≥ Cgap p3/4 , we have
(0)
sin Θ(Ûk , Uk )

 √p λ + (p p p )1/4 
1 2 3
k
≤ c0
λ2
23

(21)

with probability at least 1 − C exp(−cp). The proof of this step is closely related to the proof
for Theorem 3 in [26]. Note that
(0)

Û1

= SVDr1 (Y1 ) ,

Y1 = X1 + Z1 ,

where X1 is a fixed matrix satisfying rank(X1 ) = r1 ; Z1 ∈ Rp1 ×(p2 p3 ) ,
iid

1 ,p2 p3
∼ N (0, 1). This shares the same setting as the one in Theorem 3 in [26], if one
{(Z1 )ij }pi,j=1

sets p1 , p2 in the statement of Theorem 3 in [26] respectively as p2 p3 , p1 in our context. Thus
we can essentially follow their proof. Let U1⊥ be the orthogonal complement of U1 . Then the
Appendix Equations (1.15), (1.16) in [26] yields
2
(0)
sin Θ(Û1 , U1 )

P

≤

λ4
≥1 − C exp −c 2
λ + p2 p3


C(λ2 + p2 p3 )kU1⊥ Y1 PU > Y1 k2

!

1

λ4

(22)



and


P kU1⊥ Y1 PU > Y1 k ≥ x
1
o
n

p

+ C exp −c(λ2 + p2 p3 )
≤C exp Cp1 − c min x2 , x λ2 + p2 p3

(23)

for some uniform constant C, c > 0. Since


1/2
λ ≥ Cgap p3/4 ≥ Cgap c p1 + (p1 p2 p3 )1/4 ,
√
if we set x = C p1 , (23) further leads to

√ 
P kU1⊥ Y1 PU > Y1 k ≥ C p1 ≤ Ce−cp + C exp(−c(λ2 + p2 )) ≤ C exp(−cp).
1

(24)

Combining (24) and (22), we have proved (21) for k = 1. The proof for (21) for k = 2, 3 can
be similarly written down.
2. After spectral initialization, we assume the algorithm evolves from t = 0 to t = tmax , where

(t
)
(t
)
(t
)
tmax ≥ C log( λp ) ∨ 1 . In this step, we derive the perturbation bounds for Û1 max , Û2 max , Û3 max
3/4

under the assumptions that λ ≥ Cgap for large constant Cgap > 0 and the following inequalities all holds,
max

n

(0)

(0)

(0)

sin Θ(Û1 , U1 ) , sin Θ(Û2 , U2 ) , sin Θ(Û3 , U3 )
24

o

1
≤ ,
2

(25)

max

kZ1 · (V2 ⊗ V3 )k
√
≤ C1 pr,
kV2 kkV3 k

max

kZ3 · (V1 ⊗ V2 )k
√
≤ C1 pr,
kV1 kkV2 k

V2 ∈Rp2 ×r2
V3 ∈Rp3 ×r3

V1 ∈Rp1 ×r1
V2 ∈Rp2 ×r2

max

V3 ∈Rp3 ×r3
V1 ∈Rp1 ×r1

kZ2 · (V3 ⊗ V1 )k
√
≤ C1 pr,
kV3 kkV1 k
(26)

√
√
√
kZ1 (U2 ⊗ U3 )k ≤ C1 p1 , kZ2 (U3 ⊗ U1 )k ≤ C2 p2 , kZ3 (U1 ⊗ U2 )k ≤ C2 p3 .

(27)

Recall here that U1 , U2 , U3 are left singular subspaces for X1 , X2 , X3 , respectively. We let Lt
(t)

be the spectral sin Θ norm error for Ûk ,
Lt = max

k=1,2,3

(t)

sin Θ(Ûk , Uk ) ,

t = 0, 1, 2, . . .

(28)

Given (25), L0 ≤ 12 . Next we aim to prove that for t = 0, 1, . . .,
Lt+1 =

(t+1)
max sin Θ(Ûk
, Uk )
k=1,2,3

√
√
C1 pr
C2 p
1
≤
Lt +
≤ .
λ
λ
2
(t+1)

To show (29), we first focus on the upper bound of k sin Θ(Û1

(29)

, U1 )k when t = 0. Define

the following key components in our analysis as follows,




(t) >
(t) > Lemma 4
(t)
(t)
= M1 Y ×2 Û2
×3 Û3
= Y1 · Û2 ⊗ Û3
∈ Rp1 ×r2 r3 ,




(t)
(t) >
(t) > Lemma 4
(t)
(t)
X1 = M1 X ×2 Û2
×3 Û3
= X1 · Û3 ⊗ Û1
∈ Rp1 ×r2 r3 ,




(t)
(t) >
(t) > Lemma 4
(t)
(t)
Z1 = M1 Z ×2 Û2
×3 Û3
= Z1 · Û1 ⊗ Û2
∈ Rp1 ×r2 r3 .
(t)

Y1

By definition, the left and right singular subspaces of X1 are U1 ∈ Op1 ,r1 and U2 ⊗ U3 ∈
Op2 p3 ,r2 r3 . Then,








(t)
(t)
(t)
(t)
(t)
σr1 X1 = σr1 X1 · Û2 ⊗ Û3
= σr1 X1 · PU2 ⊗U3 · Û2 ⊗ Û3



(t)
(t)
=σr1 X1 · (U2 ⊗ U3 ) · (U2 ⊗ U3 )> · Û2 ⊗ Û3



(t)
(t)
>
≥σr1 (X1 · (U2 ⊗ U3 )) · σmin (U2 ⊗ U3 ) · Û2 ⊗ Û3

 

(t)
(t)
=σr1 (X1 ) · σmin U2> Û2 ⊗ U3> Û3




(t)
(t)
≥σr1 (X1 ) · σmin U2> Û2 · σmin U3> Û3
≥λ · 1 − L2t



(by (28) and Lemma 1 in [26]).

25

(30)

Meanwhile,
(t)

Z1



(t)
(t)
= Z1 Û2 ⊗ Û3


  (t)
(t)
= Z1 PU2 ⊗U3 + PU2⊥ ⊗U3 + PIp ⊗U3⊥ Û2 ⊗ Û3
(by Lemma 4)




(t)
(t)
(t)
(t)
+ Z1 (PU2⊥ ⊗U3 ) Û2 ⊗ Û3
≤ Z1 (PU2 ⊗U3 ) Û2 ⊗ Û3




(t)
(t)
(t)
(t)
+ Z1 (PU2⊥ ⊗U3⊥ ) Û2 ⊗ Û3
+ Z1 (PU2 ⊗U3⊥ ) Û2 ⊗ Û3


(t)
(t)
= Z1 (U2 ⊗ U3 ) (U2 ⊗ U3 )> Û2 ⊗ Û3


(t)
(t)
+ Z1 (PU2⊥ Û2 ) ⊗ (PU3 Û3 )




(t)
(t)
(t)
(t)
+ Z1 (PU2 Û2 ) ⊗ (PU3⊥ Û3 ) + Z1 (PU2⊥ Û2 ) ⊗ (PU3⊥ Û3 )
(26)

√
(t)
≤ kZ1 (U2 ⊗ U3 )k + C1 pr PU2⊥ Û2
√
(t)
+ C1 pr PU2 Û2

(t)

PU3⊥ Û3

(31)

(t)

PU3 Û3

√
(t)
+ C1 pr PU2⊥ Û2

(t)

PU3⊥ Û3

(27)

√
√
√
√
≤ C2 p1 + C1 prLt + C1 prLt + C1 prL2t
(by (27) and Lemma 1 in [26])
√
√
≤C2 p1 + 3C1 prLt
(t+1)

Since U1 and Û1

(since the spectral sin Θ norm is at most 1).
(t)

(t)

are respectively the leading r singular vectors of X1 and Y1 , by Wedin’s

sin Θ theorem [61],
sin Θ
(25) 2C2

≤

√



p1

λ
We can similarly prove that

(t+1)
Û1
, U1



√
√
C2 p1 + C1 prLt
≤
≤
(t) 
λ(1 − L2t )
σr X1
(t)

kZ1 k

√
4C1 pr
+
Lt .
λ

(32)

√
√


2C2 p2 4C1 pr
(t+1)
sin Θ Û2
, U2 ≤
+
Lt ,
λ√
λ√


2C2 p3 4C1 pr
(t+1)
sin Θ Û3
, U3 ≤
+
Lt .
λ
λ
Finally, since max{p1 , p2 , p3 } ≤ C0 p and max{r1 , r2 , r3 } ≤ C0 p1/2 , there exists a large constant Cgap > 0 such that when λ ≥ Cgap p3/4 ,
√
√
2C2 p1 4C1 pr
1
+
Lt ≤
λ
λ
2

and

√
4C1 pr
1
≤
λ
2

(33)

Then we have finished the proof for (29) for t = 0. By induction, we can sequentially prove
that (29) for all t ≥ 0.
26

At this point, (29) yields
√
√
C p 4C1 pr
Lt+1 ≤
+
Lt ,
t = 1, 2, . . . , tmax − 1
λ√
λ√ 
√ 
2C p
4C1 pr
2C p
⇒ Lt+1 −
≤
Lt −
, (since (33)),
λ
λ
λ


√
√ 
√ 
2C p
4C1 pr tmax
2C p
⇒ Ltmax −
≤
· L0 −
λ
λ
λ
!
√
√
√
√
2C p
3C p
2C p
pλ + p3/2
L0
1
≤
⇒ Ltmax ≤
+ tmax =
+ tmax C
λ
2
λ
2
λ2
λ
when tmax ≥ C log

p
λ


∨ 1 . Therefore, we have the following upper bound for spectral sin Θ

norm loss for Ûktmax = Ûk ,
k sin Θ(Ûk , Uk )k ≤ Ltmax

√
C pk
≤
.
λ

(34)

when (25), (26), (27) holds.
√
By the same calculation, we can also prove Ltmax −1 satisfies Ltmax −1 ≤ C p/λ. We prepare
this inequality for the use in the next step.
3. In this step, we develop the upper bound for X̂ − X

F

under the assumptions of (25), (26),

(27), and
Z ×1 Û1> ×2 Û2> ×3 Û3>
(t

F

√
√
√
≤ C ( p1 r1 + p2 r2 + p3 r3 ) .

(35)

)

Instead of working on kXkF and Ûk max directly, we take one step back and work on the
(t

evolution of Ûk max

−1)

(t

)

to Ûk max .

(tmax )

Recall that Û1 = Û1

(tmax )

, Û2 = Û2

(tmax )

, Û3 = Û3

; Û1⊥ , Û2⊥ , Û3⊥ are the orthogonal

complements of Û1 , Û2 , Û3 , respectively; X̂ = Y ×1 PÛ1 ×2 PÛ2 ×3 PÛ3 . In the previous step
we have also proved that
√
sin Θ(Ûk , Uk ) ≤ C

27

pk
,
λ

k = 1, 2, 3.

Then we have the following decomposition for the estimation error
X̂ − X

F

≤ X − X ×1 PÛ1 ×2 PÛ2 ×3 PÛ3

+ Z ×1 PÛ1 ×2 PÛ2 ×3 PÛ3

F

= X ×1 PÛ1⊥ + X ×1 PÛ1 ×2 PÛ2⊥ + X ×1 PÛ1 ×2 PUˆ2 ×3 PÛ3⊥

F

+ Z
≤ X

×1 Û1>

>
×1 Û1⊥

×2 Û2>

F

+ X

×3 Û3>
>
×2 Û2⊥

+ Z ×1 Û1> ×2 Û2> ×3 Û3>

F

(36)

F
>
+ X ×3 Û3⊥

F
F

F

.

To obtain the upper bound of kX̂ − Xk, we only need to analyze the four terms in (36)
separately. Recall in Step 2, we defined


(t
−1)
(t
−1)
M1 Y ×2 (Û2 max )> ×3 (Û3 max )>


(t
−1)
(t
−1)
(t
−1)
⊗ Û3 max
:= Y1 max ,
=Y1 · Û2 max
(t

X1 max

−1)

(t

, Z1 max

−1)

are defined similarly. Based on the calculation in (30) and (31), we have




 3
(t
−1)
(t
−1)
) ≥ σr (X1 ) · σmin U2> Û2 max
· σmin U3> Û3 max
≥ λ,
4
√
p1
√
√
√
√
√
(t
−1)
kZ1 max k ≤ C p1 + C prLtmax −1 ≤ C p1 + C pr ·
≤ C p1 .
λ
(t

σmin (X1 max

−1)

(tmax −1)

Since Û1 is the leading r left singular vectors of Y1

(t

= X1 max

−1)

(t

+ Z1 max

−1)

, Lemma 6

implies


(t
−1)
(t
−1)
PÛ1⊥ M1 X ×2 (Û2 max )> ×3 (Û3 max )>
=

(t
−1)
PÛ1⊥ X1 max

F

√

F

≤ C p1 r1 .

As a result,
X ×1 PÛ1⊥ = PÛ1⊥ · X1 · (PU2 ⊗ PU3 ) = PÛ1⊥ · X1 · (U2 ⊗ U3 )
F
F
F


(t
−1)
(t
−1)
(t
−1)
(t
−1)
−1
−1
≤ PÛ1⊥ X1 Û2 max
⊗ Û3 max
· σmin
(U2> Û2 max ) · σmin
(U3> Û3 max )
F

(37)

1
1
√
√
p
≤C p1 r1 p
≤ C p1 r1 .
2
2
1 − (1/2)
1 − (1/2)
Similarly, we can show
X̂ ×2 PÛ2⊥

F

√
≤ C p2 r2 ,
28

X̂ ×3 PÛ3⊥

F

√
≤ C p3 r3 .

(38)

Now combining (36), (35), (37), and (38), we have
X̂ − X

F

√
√
√
≤ C ( p1 r1 + p2 r2 + p3 r3 )

(39)

for some constant C > 0.
4. We finalize the proof for Theorem 1 in this step. By Lemma 5, we know (26), (27), and (35)
hold with probability at least 1 − C exp(−cp). By the result in Step 1, we know (25) holds
with probability at least 1 − C exp(−cp). Let Q = {(26), (27), (35), (25) all hold}, then
P (Q) ≥ 1 − C exp(−cp).

(40)

√
By Steps 2 and 3, one has k sin Θ(Ûk , Uk )k ≤ C pk /λ, k = 1, 2, 3, and
√
√
√
kX̂ − XkF ≤ C ( p1 r1 + p2 r2 + p3 r3 ) under Q.
It remains to consider situation under Qc . By definition, X̂ is a projection of Y, so
kX̂kF ≤ kYkF ≤ kXkF + kZkF .
Then we have the following rough upper bound for the 4-th moment of recovery error,


EkX̂ − Xk4F ≤ C EkX̂k4F + kXk4F ≤ CkXk4F + CEkZk4F
≤C exp(c0 p) + CE χ2p1 p2 p3

2

≤ C exp(c0 p) + Cp6 .

The the following upper bound holds for the Frobenius norm risk of X̂,
EkX̂ − Xk2F = EkX̂ − Xk2F 1Q + EkX̂ − Xk2F 1Qc
q
√
√
√
=C ( p1 r1 + p2 r2 + p3 r3 ) + EkX̂ − Xk4F · EQc
(40)

√
√
√
≤ C ( p1 r1 + p2 r2 + p3 r3 ) + C exp ((c0 − c)p/2) + Cp3 exp(−cp/2).

Thus, one can select c0 < c to ensure that
EkX̂ − Xk2F ≤ C (p1 r1 + p2 r2 + p3 r3 ) .
Additionally, since σrk (Mk (X)) ≥ λ, we have kXkF2 = kMk (X)k2F ≥ rk λ2 for k = 1, 2, 3,
which implies kXk2F ≥ max{r1 , r2 , r3 }λ = λr, then
E

kX̂ − Xk2F
p1 + p2 + p3
≤
.
2
λ
kXkF
29

Moreover, by definition, k sin Θ(Ûk , Uk )k ≤ 1. Thus we have the following upper bound for
the spectral sin Θ risk for Ûk ,
Ek sin Θ(Ûk , Uk )k ≤ Ek sin Θ(Ûk , Uk )k1Q + Ek sin Θ(Ûk , Uk )k1Qc
√
√
q
(40)
pk
pk p
4
=C
+ Ek sin Θ(Ûk , Uk )k · E1Qc ≤ C
+ C exp(−cp).
λ
λ
By definition of λ, we know λ = σrk (Mk (X)) ≤

kXkF
√
rk

≤

C exp(c0 p)
√
,
rk

so one can select small

constant c0 > 0 to ensure that
√

√
p
pk
pk rk
≥
≥ c exp(−cp),
λ
C exp(c0 p)

which implies Ek sin Θ(Ûk , Uk )k ≤ C pλk . Finally, we can derive the general Schatten q-sin Θnorm risk via Hölder’s inequality,
√
−1/q
Erk Ek sin Θ(Ûk , Uk )kq

≤ Ek sin Θ(Ûk , Uk )k ≤ C

pk
.
λ

Summarizing from Steps 1-4, we have finished the proof of Theorem 1.

7.2



Proof of Theorem 4

We particularly show that it suffices to consider sparse tensor models and we set σ = 1 for
brevity. A tensor X ∈ Rp1 ×p2 ×p3 is sparse with respect to parameters S(X) = (s1 , s2 , s3 ) if there
exists Sk (X) ⊂ [pk ] := {1, 2, . . . , pk }, k = 1, 2, 3 such that
Xijk = 0,

∀(i, j, k) ∈ [p1 ] × [p2 ] × [p3 ] \ S1 (X) × S2 (X) × S3 (X)

with |Sk (X)| ≤ sk , k = 1, 2, 3. It means that the nonzero entries of X are constrained in the block
S1 (X) × S2 (X) × S3 (X). Define the subset M(p, k, r, λ) ⊂ Fp,r (λ) for integer k = bp(1−τ )/2 c as
follows,
n
o
M(p, k, r, λ) := X ∈ Fp,r (λ) : S(X) ≤ (20k, 20k, 20k) ,
containing sparse tensors in Fp,r (λ). Consider two disjoint subsets of M(p, k, r, λ):
n
o
M0 (p, k, r, λ) := X ∈ M(p, k, r, λ), S1 (X) ∪ S2 (X) ∪ S3 (X) ⊂ [p/2] ,
and
n
o
M1 (p, k, r, λ) := X ∈ M(p, k, r, λ), S1 (X) ∪ S2 (X) ∪ S3 (X) ⊂ [p] \ [p/2]
30

where matrices in M0 (p, k, r, λ) and M1 (p, k, r, λ) are supported on disjoint blocks, so are their
singular vectors. Given the observation:
Y = X + Z ∈ Rp1 ×p2 ×p3 ,
the following testing problem is studied:
H0 : X ∈ M0 (p, k, r, λ)

V.S. H1 : X ∈ M1 (p, k, r, λ).

(41)

A test is then defined as φ(·) : Rp1 ×p2 ×p3 → {0, 1} whose risk is given as
Rp,r,λ (φ) =


PX φ(Y) = 1 +

sup
X∈M0 (p,k,r,λ)


PX φ(Y) = 0 ,

sup

X∈M1 (p,k,r,λ)

the worst case of Type-I+II error.
Lemma 1. Suppose Hypothesis H(τ ) for some τ ∈ (0, 1). Let {φp } be any sequence of polynomial 3(1−τ )/4 
time tests of (41). There exists an absolute constant c0 > 0 such that if λ ≤ c0 p √log p , then
as long as min{r1 , r2 , r3 } ≥ 1,
1
lim inf Rp,r,λ (φp ) ≥ .
p→∞
2
Now we move back to the proof for Theorem 4. Suppose that, on the contradiction, for any
(p)

k = 1, 2, 3, there exists a sub-sequence (Ûk ) such that
lim


(p)
E sin Θ Ûk , Uk (X) = 0,

sup

p→∞ X∈F

p,r (λ)

which implies that
lim



sup

p→∞ X∈F

P

p,r (λ)

(p)

Uk Uk> − Ûk

(p) >

≤

Ûk

1
= 1.
3

Define a sequence of tests φp : Rp1 ×p2 ×p3 7→ {0, 1} as

φp (Y) :=

(p) 

where Ûk



0,

if


1,

otherwise,

(p) 

Ûk

[1:p/2,:]

(p) >
[1:p/2,:]

Ûk

≥ 23 ,

(p)

[1:p/2,:]

denote the first p/2 rows of Ûk . Clearly,

Rp,r,λ (φp ) ≤

sup


PX

X∈M0 (p,k,r,λ)

(p)

Uk Uk> − Ûk

31

(p) >

Ûk

>

1
3

(42)

+



sup

PX

X∈M1 (p,k,r,λ)

which implies limp→∞ Rp,r,λ φp



(p)

Uk Uk> − Ûk

(p) >

Ûk

≥

2
,
3

= 0, contradicting Lemma 1. Now, we prove claim (18).

Suppose that, on the contradiction, there exists a sub-sequence (X̂(p) ) such that
lim

sup

p→∞ X∈F

E

p,r (λ)

kX̂(p) − Xk2F
= 0,
kXk2F

which implies


1
P kX̂(p) − XkF ≤ kXkF = 1.
p→∞ X∈F (λ)
3
p,r
lim

(43)

sup

Define a sequence of test φp : Rp1 ×p2 ×p3 7→ {0, 1} as

φp (Y) :=



0,

if


1,

otherwise,

X̂(p)


[V1 ,V1 ,V1 ] F

≥

X̂(p)


[V2 ,V2 ,V2 ] F

,

where V1 = [p/2], V2 = [p] \ V1 and X[V1 ,V1 ,V1 ] denotes the sub-tensor on the block V1 × V1 × V1 .
Under H0 , if φp (Y) = 1, then
kX̂(p) − Xk2F ≥
≥

X̂(p) − X

X̂(p) − X




2
[V1 ,V1 ,V1 ] F

2
[V1 ,V1 ,V1 ] F

X̂(p)

+



+

X̂(p)



2
[V1 ,V1 ,V1 ] F

2
[V2 ,V2 ,V2 ] F

1
≥ kXk2F .
2

Clearly,
√

2
kXkF
Rp,r,λ (φp ) =
sup
P kX̂ − XkF ≥
2
X∈M0 (p,k,r,λ)
√


2
kXkF ,
+
sup
P kX̂(p) − XkF ≥
2
X∈M1 (p,k,r,λ)


(p)

which implies that Rp,r,λ (φp ) → 0 as p → ∞ based on (43), which contradicts Lemma 1.

7.3



Proof of Theorems 3

Without loss of generality we can assume σ = 1 throughout the proof. First, we construct
the core tensor S̃ ∈ Rr1 ×r2 ×r3 with i.i.d. standard Gaussian entries, then according to random
matrix theory (c.f. Corollary 5.35 in [62]), with probability at least 1 − 6e−x , we have
√
√
√
√
rk+1 rk+2 − rk − x ≤ σmin (Mk (S̃)) ≤ σmax (Mk (S̃)) ≤ rk+1 rk+2 + rk + x,
32

for k = 1, 2, 3. Plug in x = 1.8, by simple calculation, we can see there is a positive probability
that
√

rk+1 rk+2 −

√
√
√
rk − 1.8 ≤ σmin (Mk (S̃)) ≤ σmax (Mk (S̃)) ≤ rk+1 rk+2 + rk + 1.8.

(44)

Note that
r1 r2 ≥ 4r3 , r2 r3 ≥ 4r1 , r3 r1 ≥ 4r2 ,

⇒ r1 r2 r3 ≥ 4 max {rk }2
1≤k≤3

rk+1
rk+2
≥ 4, ⇒ rk ≥ 4
max{rk } max{rk }
4 max{rk }2
⇒ rk+1 rk+2 ≥
≥ 4rk ≥ 16,
rk
⇒ rk

we know
√
√
rk+1 rk+2 − rk − 1.8
r
rk+1 rk2
1.8 √
√
√
−
≥ rk+1 rk+2 −
rk+1 rk+2 = 0.05 rk+1 rk+2 ;
4
4
√
√
rk+1 rk+2 − rk + 1.8
r
1.8 √
rk+1 rk2
√
√
+
rk+1 rk+2 = 1.95 rk+1 rk+2 .
≤ rk+1 rk+2 +
4
4
√
By previous arguments, there exists S̃ ∈ Rr1 ×r2 ×r3 such that c rk+1 rk+2 ≤ σmin (Mk (S̃)) ≤
√
λ
C rk+1 rk+2 for k = 1, 2, 3. Now, we construct the scaled core tensor S = S̃ min
.
σ
(M (S̃))
k=1,2,3

min

k

Given r ≤ r1 , r2 , r3 ≤ C0 r, we know S ∈ Rr1 ×r2 ×r3 satisfies the following property
λ ≤ σmin (Mk (S)) ≤ σmax (Mk (S)) ≤ Cλ,

Proof of the first claim.

k = 1, 2, 3.

(45)

It suffices to consider k = 1. We construct a large subset of Op1 ,r1

whose elements are well separated in Schatten q-norms for all 1 ≤ q ≤ +∞. To this end, we
need some preliminary facts about the packing number in Grassmann manifold Gp,r , which is
the set of all r-dimensional subspaces of Rp . Given such a subspace L ⊂ Rp with dim(L) = r,
let UL ∈ O(p, r) denote the orthonormal basis of L. Denote Bp,r := {UL , L ∈ Gp,r } which is
actually a subset of Op,r and will be equipped with Schatten q-norm distances for all q ∈ [1, +∞]:
dq (UL1 , UL2 ) := kUL1 UL>1 − UL2 UL>2 kq . Recall that the ε-packing number of a metric space (T, d)
is defined as
n
o
D(T, d, ε) := max n : there are t1 , . . . , tn ∈ T, such that min d(ti , tj ) > ε .
i6=j

33

The following lemma can be found in Lemma 5 in [63] which controls the packing numbers of
Bp,r with respect to Schatten distances dq .
Lemma 2. For all integers 1 ≤ r ≤ p such that r ≤ p − r, and all 1 ≤ q ≤ +∞, the following
bound holds
 c r(p−r)
ε

  C r(p−r)
≤ D Bp,r , dq , εr1/q ≤
ε

with absolute constants c, C > 0.
We are in position to construct a well-separated subset of Op1 ,r1 . According to Lemma 2 by

choosing ε = 2c , there exists a subset Vp1 −r1 ,r1 ⊂ Op1 −r1 ,r1 with Card Vp1 −r1 ,r1 ≥ 2r1 (p1 −2r1 )
such that for each V1 6= V2 ∈ Vp1 −r1 ,r1 ,
c 1/q
kV1 V1> − V2 V2> kq ≥ r1 .
2
Now, fix a δ > 0 whose value is to be determined later. For every V ∈ Vp1 −r1 ,r1 , define Ṽ ∈ Op1 ,r1
as follows

 √

1 − δIr1
.
Ṽ = 
√
δV

It is easy to check that Ṽ ∈ Op1 ,r1 as long as V ∈ Op1 −r1 ,r1 . We conclude with a subset Vp1 ,r1 ⊂
Op1 ,r1 with Card(Vp1 ,r1 ) = Card(Vp1 −r1 ,r1 ) ≥ 2r1 (p1 −2r1 ) . Moreover, for Ṽ1 6= Ṽ2 ∈ Vp1 ,r1 ,
kṼ1 Ṽ1> − Ṽ2 Ṽ2> kq ≥

p
cp
1/q
δ(1 − δ)kV1 − V2 kq ≥
δ(1 − δ)r1 .
2

Meanwhile,
kṼ1 − Ṽ2 kF ≤

√

δkV1 − V2 kF ≤

p
2δr1 .

Then we construct a series of fixed signal tensors: Xi = S ×1 Ṽi ×2 U2 ×3 U3 , i = 1, . . . , m,
where U2 ∈ Op2 ,r2 .U3 ∈ Op3 ,r3 are any fixed orthonormal columns, Ṽi ∈ Vp1 ,r1 ⊂ Op1 ,r1 and m =
2 (M (X )) ≥ λ for k = 1, 2, 3, so that {X }m ⊆ F
2r1 (p1 −2r1 ) . By such construction, σmin
i
i i=1
p,r (λ).
k

We further let Yi = Xi + Zi , where Zi are i.i.d. standard normal distributed tensors, which
implies Yi ∼ N (Xi , Ip1 ×p2 ×p3 ). Then the Kullback-Leibler divergence between the distribution
Yi and Yj is
DKL (Yi ||Yj ) =

1
1
kXi − Xj k2F =
S ×1 (Ṽi − Ṽj ) ×2 U2 ×3 U3
2
2

1
= k(Ṽi − Ṽj ) · M1 (S) · (U2 ⊗ U3 )> k2F ≤ Cλ2 Ṽi − Ṽj
2
34

2
F

2

2
F

≤ 2Cλ δr1 .

Then the generalized Fano’s lemma yields the following lower bound


cp
Cλ2 δr1 + log 2
1/q
sup E Û1 Û1> − U1 U1> ≥
inf
δ(1 − δ)r1
1−
.
2
r
(p
−
2r
)
log
2
q
1
1
1
Û1 U1 ∈{Vi }m
i=1

(46)

1)
By setting δ = c1 (p1 −2r
for a small but absolute constant c1 > 0, we obtain
λ2



Cλ2 δr1 + log 2
(p1 − 2r1 )1/2 1/q
cp
1/q
δ(1 − δ)r1
1−
≥ c0
r1
2
r1 (p1 − 2r1 ) log 2
λ
for an absolute constant c0 > 0. Then, if p1 ≥ 3r1 ,
sup

inf

Û1 U1 ∈{Vi }m
i=1
1/q

where r1

E

Û1 Û1>

−

U1 U1>

q

≥ c0

 √p
λ

1 1/q
r1

1/q

∧ r1



is a trivial term. The first claim in Theorem 3 it thus obtained by viewing the

equivalence between the Schatten q-norms and sin Θ Schatten q-norms, see Lemma 3 in the
Appendix.

Proof of second and third claims.

To prove the minimax lower bounds in estimating X, we

need a different construction scheme. Specifically, we consider the metric space (Op1 ,r1 , k sin Θ(·, ·)k2 ),
fix an V0 ∈ Op1 ,r1 , and consider the following ball of radius ε > 0 and center V0 :

B(V0 , ε) = V 0 ∈ Op1 ,r1 : k sin Θ(V 0 , V )k2 ≤ ε .
By Lemma 1 in [10], for 0 < α < 1 and 0 < ε ≤ 1, there exists V10 , . . . , Vm0 ⊆ B(V0 , ε) such that
m≥

 c r1 (p1 −r1 )
0

α

,

min

1≤i<j≤m

sin Θ(Vi0 , Vj0 )

2

≥ αε.

By the property of sin Θ distance (Lemma 1 in [26]), we can find a rotation matrix Oi ∈ Or1
such that
kV0 − Vi0 Oi kF ≤

√

2k sin Θ(V0 , Vi0 )k2 ≤

√
2ε.

We denote Vi = Vi0 Oi , then
kVi − V0 kF ≤

√

2ε,

k sin Θ(Vi , Vj )k2 ≥ αε,

35

1 ≤ i < j ≤ m.

(47)

Construct Xi = S ×1 Vi ×2 U2 × U3 for i = 1, . . . , m in a similar fashion as above. Then the class
of low-rank tensors satisfy the following properties,
kXi − Xj k2F = kS ×1 (Vi − Vj ) ×2 U2 ×3 U3 k2F
1
1
= k(Vi − Vj ) · M1 (S) · (U2 ⊗ U3 )> k2F ≥ σr1 (M1 (S)) kVi − Vj k2F
2
2
2
λ
≥
min kVi − Vj Ok2F (by Lemma 1 in [26])
2 O∈Or
λ2
≥ ksin Θ (Vi , Vj )k22 ≥ α2 ε2 λ2 , 1 ≤ i < j ≤ m.
2

(48)

√
min kXi kF = kS ×1 Vi ×2 U2 ×3 U3 kF = kSkF ≥ kM1 (S)kF ≥ λ r1 .

(49)

1≤i≤m

Moreover, under the same model Yi = Xi + Zi as above, the KL-divergence between the
distributions of Yi and Yj is
DKL (Yi ||Yj ) =

1
1
kXi − Xj k2F = kS ×1 (Vi − Vj ) ×2 U2 ×3 U3 k2F
2
2

1
= k(Vi − Vj ) · M1 (S) · (U2 ⊗ U3 )> k2F ≤ Cλ2 kVi − Vj k2F ≤ Cλ2 ε2 .
2
Then the generalized Fano’s lemma yields the following lower bound


2
Cλ2 ε2 + log 2
inf
sup E X̂ − X ≥ λαε 1 −
.
r
(p
−
r
)
log(c
/α)
F
1
1
1
0
X̂ X∈{Xi }m
i=1
By setting ε =

q

r1 (p1 −r1 )
2Cλ2

∧

√


αε 1 −

(50)

2r1 , α = (c0 ∧ 1)/8, we have

Cλ2 ε2 + log 2
r1 (p1 − r1 ) log(c0 /α)

√


≥ c1

r1 p1 √
∧ r1
λ



for some small constant c1 > 0. Moreover,
inf

E X̂ − X

sup

X̂ X∈Fp,r (λ)

≥ inf

sup

X̂ X∈{X1 ,...,Xm }

2
F

≥ inf

E X̂ − X

sup

√
X̂ X∈Fp,r (λ∨ p1 )

2 (50)
F

2

≥ c1 (λ ∨ p1 )

E X̂ − X



2
F

r1 p1
∧ r1
2
λ ∨ p1

(51)


≥ c1 p1 r1 .

kX̂ − Xk2F
kX̂ − Xk2F
≥
inf
sup
E
kXk2F
kXk2F
X̂ X∈Fp,r (λ)
X̂ X∈Fp,r (λ∨p1 )

p

(49)(50) c λ2  r p
kX̂ − Xk2F
1
1 1
1
≥ inf
sup
E
≥
≥
c
∧
r
∧
1
.
1
1
2
λ2 r1 λ2
λ2
X̂ X∈{X1 ,...,Xm } max1≤i≤m kXkF
inf

sup

E

(52)

Finally, we apply the same argument of (50), (51), and (52) to U2 , U3 , then we can obtain
(13) and (14).


36

Acknowledgment
The authors would like to thank Ming Yuan and Zongming Ma for helpful discussions. The
authors would also like to thank Genevera Allen and Nathaniel Helwig for pointing out several
key references. The authors would also like to thank editors and anomalous referees for your
helpful comments and suggestions on improving the paper.

References
[1] H. Shen and J. Z. Huang, “Sparse principal component analysis via regularized low rank
matrix approximation,” Journal of multivariate analysis, vol. 99, no. 6, pp. 1015–1034,
2008.
[2] M. Lee, H. Shen, J. Z. Huang, and J. Marron, “Biclustering via sparse singular value
decomposition,” Biometrics, vol. 66, no. 4, pp. 1087–1095, 2010.
[3] D. Yang, Z. Ma, and A. Buja, “A sparse singular value decomposition method for highdimensional data,” Journal of Computational and Graphical Statistics, vol. 23, no. 4,
pp. 923–942, 2014.
[4] D. Yang, Z. Ma, and A. Buja, “Rate optimal denoising of simultaneously sparse and low
rank matrices,” Journal of Machine Learning Research, vol. 17, no. 92, pp. 1–27, 2016.
[5] E. J. Candes, C. A. Sing-Long, and J. D. Trzasko, “Unbiased risk estimates for singular
value thresholding and spectral estimators,” IEEE transactions on signal processing, vol. 61,
no. 19, pp. 4643–4657, 2013.
[6] A. A. Shabalin and A. B. Nobel, “Reconstruction of a low-rank matrix in the presence of
gaussian noise,” Journal of Multivariate Analysis, vol. 118, pp. 67–76, 2013.
[7] D. Donoho, M. Gavish, et al., “Minimax risk of matrix denoising by singular value thresholding,” The Annals of Statistics, vol. 42, no. 6, pp. 2413–2440, 2014.
√
[8] M. Gavish and D. L. Donoho, “The optimal hard threshold for singular values is 4/ 3,”
IEEE Transactions on Information Theory, vol. 60, no. 8, pp. 5040–5053, 2014.
37

[9] H. Zou, T. Hastie, and R. Tibshirani, “Sparse principal component analysis,” Journal of
computational and graphical statistics, vol. 15, no. 2, pp. 265–286, 2006.
[10] T. T. Cai, Z. Ma, Y. Wu, et al., “Sparse pca: Optimal rates and adaptive estimation,” The
Annals of Statistics, vol. 41, no. 6, pp. 3074–3110, 2013.
[11] A. Birnbaum, I. M. Johnstone, B. Nadler, and D. Paul, “Minimax bounds for sparse pca
with noisy high-dimensional data,” Annals of statistics, vol. 41, no. 3, p. 1055, 2013.
[12] E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis?,” Journal
of the ACM (JACM), vol. 58, no. 3, p. 11, 2011.
[13] H. Zhou, L. Li, and H. Zhu, “Tensor regression with applications in neuroimaging data
analysis,” Journal of the American Statistical Association, vol. 108, no. 502, pp. 540–552,
2013.
[14] A. Zhang, “Cross: Efficient low-rank tensor completion,” arXiv preprint arXiv:1611.01129,
2016.
[15] A. Karatzoglou, X. Amatriain, L. Baltrunas, and N. Oliver, “Multiverse recommendation:
n-dimensional tensor factorization for context-aware collaborative filtering,” in Proceedings
of the fourth ACM conference on Recommender systems, pp. 79–86, ACM, 2010.
[16] S. Rendle and L. Schmidt-Thieme, “Pairwise interaction tensor factorization for personalized tag recommendation,” in Proceedings of the third ACM international conference on
Web search and data mining, pp. 81–90, ACM, 2010.
[17] H. Lu, K. N. Plataniotis, and A. N. Venetsanopoulos, “Mpca: Multilinear principal component analysis of tensor objects,” IEEE Transactions on Neural Networks, vol. 19, no. 1,
pp. 18–39, 2008.
[18] J. Liu, P. Musialski, P. Wonka, and J. Ye, “Tensor completion for estimating missing values
in visual data,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 35,
no. 1, pp. 208–220, 2013.

38

[19] A. Anandkumar, R. Ge, D. J. Hsu, S. M. Kakade, and M. Telgarsky, “Tensor decompositions
for learning latent variable models.,” Journal of Machine Learning Research, vol. 15, no. 1,
pp. 2773–2832, 2014.
[20] A. Anandkumar, R. Ge, D. Hsu, and S. M. Kakade, “A tensor approach to learning mixed
membership community models,” J. Mach. Learn. Res., vol. 15, pp. 2239–2312, Jan. 2014.
[21] N. Li and B. Li, “Tensor completion for on-board compression of hyperspectral images,” in
2010 IEEE International Conference on Image Processing, pp. 517–520, IEEE, 2010.
[22] T. Liu, M. Yuan, and H. Zhao, “Characterizing spatiotemporal transcriptome of human
brain via low rank tensor decomposition,” arXiv preprint arXiv:1702.07449, 2017.
[23] E. Richard and A. Montanari, “A statistical model for tensor pca,” in Advances in Neural
Information Processing Systems, pp. 2897–2905, 2014.
[24] S. B. Hopkins, J. Shi, and D. Steurer, “Tensor principal component analysis via sum-ofsquare proofs.,” in COLT, pp. 956–1006, 2015.
[25] A. Perry, A. S. Wein, and A. S. Bandeira, “Statistical limits of spiked tensor models,” arXiv
preprint arXiv:1612.07728, 2016.
[26] T. T. Cai and A. Zhang, “Rate-optimal perturbation bounds for singular subspaces with
applications to high-dimensional statistics,” The Annals of Statistics, vol. to appear, 2017.
[27] C. J. Hillar and L.-H. Lim, “Most tensor problems are np-hard,” Journal of the ACM
(JACM), vol. 60, no. 6, p. 45, 2013.
[28] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM journal on Matrix Analysis and Applications, vol. 21, no. 4, pp. 1253–1278,
2000.
[29] L. De Lathauwer, B. De Moor, and J. 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, vol. 21, no. 4, pp. 1324–1342, 2000.

39

[30] M. A. O. Vasilescu and D. Terzopoulos, “Multilinear subspace analysis of image ensembles,”
in Computer Vision and Pattern Recognition, 2003. Proceedings. 2003 IEEE Computer
Society Conference on, vol. 2, pp. II–93, IEEE, 2003.
[31] B. N. Sheehan and Y. 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, pp. 355–365, SIAM, 2007.
[32] R. Costantini, L. Sbaiz, and S. Susstrunk, “Higher order svd analysis for dynamic texture
synthesis,” IEEE Transactions on Image Processing, vol. 17, no. 1, pp. 42–52, 2008.
[33] M. Haardt, F. Roemer, and G. Del Galdo, “Higher-order svd-based subspace estimation to
improve the parameter estimation accuracy in multidimensional harmonic retrieval problems,” IEEE Transactions on Signal Processing, vol. 56, no. 7, pp. 3198–3213, 2008.
[34] Y. Liu, F. Shang, W. Fan, J. Cheng, and H. Cheng, “Generalized higher-order orthogonal
iteration for tensor decomposition and completion,” in Advances in Neural Information
Processing Systems, pp. 1763–1771, 2014.
[35] Q. Zheng and R. Tomioka, “Interpolating convex and non-convex tensor decompositions via
the subspace norm,” in Advances in Neural Information Processing Systems, pp. 3106–3113,
2015.
[36] A. Anandkumar, Y. Deng, R. Ge, and H. Mobahi, “Homotopy analysis for tensor pca,”
arXiv preprint arXiv:1610.09322, 2016.
[37] A. Anandkumar, R. Ge, and M. Janzamin, “Guaranteed non-orthogonal tensor decomposition via alternating rank-1 updates,” arXiv preprint arXiv:1402.5180, 2014.
[38] W. W. Sun, J. Lu, H. Liu, and G. Cheng, “Provable sparse tensor decomposition,” Journal
of Royal Statistical Association, 2015.
[39] G. I. Allen, “Regularized tensor factorizations and higher-order principal components analysis,” arXiv preprint arXiv:1202.2476, 2012.
[40] G. Allen, “Sparse higher-order principal components analysis.,” in AISTATS, vol. 15, 2012.

40

[41] Y. Qi, P. Comon, and L.-H. Lim, “Uniqueness of nonnegative tensor approximations,” IEEE
Transactions on Information Theory, vol. 62, no. 4, pp. 2170–2183, 2016.
[42] C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, and S. Yan, “Tensor robust principal component
analysis: Exact recovery of corrupted low-rank tensors via convex optimization,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 5249–5257,
2016.
[43] T. Lesieur, L. Miolane, M. Lelarge, F. Krzakala, and L. Zdeborová, “Statistical and computational phase transitions in spiked tensor estimation,” in Information Theory (ISIT),
2017 IEEE International Symposium on, pp. 511–515, IEEE, 2017.
[44] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM review,
vol. 51, no. 3, pp. 455–500, 2009.
[45] Q. Berthet and P. Rigollet, “Computational lower bounds for sparse pca,” arXiv preprint
arXiv:1304.0828, 2013.
[46] Q. Berthet, P. Rigollet, et al., “Optimal detection of sparse principal components in high
dimension,” The Annals of Statistics, vol. 41, no. 4, pp. 1780–1815, 2013.
[47] T. Wang, Q. Berthet, and R. J. Samworth, “Statistical and computational trade-offs in
estimation of sparse principal components,” arXiv preprint arXiv:1408.5369, 2014.
[48] C. Gao, Z. Ma, and H. H. Zhou, “Sparse cca: Adaptive estimation and computational
barriers,” arXiv preprint arXiv:1409.8565, 2014.
[49] Z. Ma and Y. Wu, “Computational barriers in minimax submatrix detection,” The Annals
of Statistics, vol. 43, no. 3, pp. 1089–1116, 2015.
[50] T. T. Cai, T. Liang, and A. Rakhlin, “Computational and statistical boundaries for submatrix localization in a large noisy matrix,” arXiv preprint arXiv:1502.01988, 2015.
[51] Y. Chen and J. Xu, “Statistical-computational tradeoffs in planted problems and submatrix localization with a growing number of clusters and submatrices,” arXiv preprint
arXiv:1402.1267, 2014.

41

[52] B. Barak and A. Moitra, “Noisy tensor completion via the sum-of-squares hierarchy,” in
29th Annual Conference on Learning Theory, pp. 417–445, 2016.
[53] B. E. Hajek, Y. Wu, and J. Xu, “Computational lower bounds for community detection on
random graphs.,” in COLT, pp. 899–928, 2015.
[54] B. Bollobás and P. Erdös, “Cliques in random graphs,” in Mathematical Proceedings of the
Cambridge Philosophical Society, vol. 80, pp. 419–427, Cambridge Univ Press, 1976.
[55] V. Feldman, E. Grigorescu, L. Reyzin, S. Vempala, and Y. Xiao, “Statistical algorithms
and a lower bound for detecting planted cliques,” in Proceedings of the forty-fifth annual
ACM symposium on Theory of computing, pp. 655–664, ACM, 2013.
[56] N. Alon, M. Krivelevich, and B. Sudakov, “Finding a large hidden clique in a random
graph,” Random Structures and Algorithms, vol. 13, no. 3-4, pp. 457–466, 1998.
[57] B. P. Ames and S. A. Vavasis, “Nuclear norm minimization for the planted clique and
biclique problems,” Mathematical programming, vol. 129, no. 1, pp. 69–89, 2011.
[58] L. Kuvcera, “Expected complexity of graph partitioning problems,” Discrete Applied Mathematics, vol. 57, no. 2, pp. 193–212, 1995.
[59] M. Jerrum, “Large cliques elude the metropolis process,” Random Structures & Algorithms,
vol. 3, no. 4, pp. 347–359, 1992.
[60] U. Feige and R. Krauthgamer, “The probable value of the lovász–schrijver relaxations for
maximum independent set,” SIAM Journal on Computing, vol. 32, no. 2, pp. 345–370, 2003.
[61] P.-A. Wedin, “Perturbation bounds in connection with singular value decomposition,” BIT
Numerical Mathematics, vol. 12, no. 1, pp. 99–111, 1972.
[62] R. Vershynin, “Introduction to the non-asymptotic analysis of random matrices,” arXiv
preprint arXiv:1011.3027, 2010.
[63] V. Koltchinskii and D. Xia, “Optimal estimation of low rank density matrices,” Journal of
Machine Learning Research, vol. 16, pp. 1757–1792, 2015.

42

[64] D. Hush and C. Scovel, “Concentration of the hypergeometric distribution,” Statistics &
probability letters, vol. 75, no. 2, pp. 127–132, 2005.
[65] B. Laurent and P. Massart, “Adaptive estimation of a quadratic functional by model selection,” Annals of Statistics, pp. 1302–1338, 2000.

43

Supplement to “Tensor SVD: Statistical and Computational
Limits” 1
Anru Zhang

and

Dong Xia

University of Wisconsin-Madison
Abstract
In this Supplement, we provide additional proofs for the main results and technical lemmas.

A

Additional Proofs

A.1

Proof of Theorem 2

We only need to prove upper bounds for Ûk• and X̂• as the ones for Ûkmle and X̂mle immediately
follow. The proof of this theorem is fairly complicated. For convenience, we assume σ = 1,
denote r = max{r1 , r2 , r3 }, then r ≤ C0 p1/2 according to the assumption. For any orthogonal
columns, e.g. Uk ∈ Opk ,rk , we note Uk⊥ ∈ Opk ,pk −rk as the orthogonal complement of Uk .
Let A ⊗ B be the Kronecker product between matrices A and B, vec(·) be vectorization of
matrices and tensors. Similarly as the proof of Theorems 1 and 2 in [29], for any Vk ∈ Opk ,rk ,
k = 1, 2, 3,
min kY − Ŝ ×1 V1 ×2 V2 ×3 V3 k2F = min kvec(Y) − V1 ⊗ V2 ⊗ V3 vec(Ŝ)k22
Ŝ

Ŝ

= min kP(V1 ⊗V2 ⊗V3 )⊥ vec(Y) + PV1 ⊗V2 ⊗V3 vec(Y) − V1 ⊗ V2 ⊗ V3 vec(Ŝ)k22
Ŝ

=kP(V1 ⊗V2 ⊗V3 )⊥ vec(Y)k22 + min kP(V1 ⊗V2 ⊗V3 ) vec(Y) − V1 ⊗ V2 ⊗ V3 vec(Ŝ)k22
Ŝ

=kP(V1 ⊗V2 ⊗V3 )⊥ vec(Y)k22 = kvec(Y)k22 − kP(V1 ⊗V2 ⊗V3 ) vec(Y)k22
=kYk2F − kY ×1 V1> ×2 V2> ×3 V3> k2F
1

Anru Zhang is Assistant Professor, Department of Statistics, University of Wisconsin-Madison, Madison,

WI 53706, E-mail: anruzhang@stat.wisc.edu; Dong Xia is Visiting Assistant Professor, Department of Statistics,
University of Wisconsin-Madison, Madison, WI 53706, E-mail: dongxia@stat.wisc.edu.

1

where the inequality holds if and only if Ŝ = Y ×1 V1> ×2 V2> ×3 V3> . Therefore, we must have
kY ×1 (Û1• )> ×2 (Û2• )> ×3 (Û3• )> k2F ≥ kY ×1 U1> ×2 U2> ×3 U3> k2F .
For convenience, we simply let Ûk = Ûk• for k = 1, 2, 3 and X̂ = X̂• throughout the proof of this
theorem. Without loss of generality, we also assume that


Irk
 , k = 1, 2, 3.
Uk = 
0(pk −rk )×rk

(53)

Such assumption will simplify our notation and make the proof easier to understand. This
theorem will be shown by steps.
1. In this first step, we establish some basic probability bounds which will be used in the latter
steps. We first let
X̃ = Y[1:r1 ,1:r2 ,1:r3 ] = S + Z[1:r1 ,1:r2 ,1:r3 ] ∈ Rr1 ×r2 ×r3 .

(54)

Then we first have
X̃

F

= Y ×1 U1> ×2 U2> ×3 U3>

F

≤ Y ×1 Û1> ×2 Û2> ×3 Û3>

F

.

(55)

Next, we also note that M1 (X̃) = M1 (S) + M1 (Z[1:r1 ,1:r2 ,1:r3 ] ) ∈ Rr1 ×(r2 r3 ) , i.e., the fixed
matrix M1 (S) plus an i.i.d. Gaussian matrix. Meanwhile, σmin (M1 (S)) = σr1 (M1 (X)) ≥ λ.
Now, by Lemma 4 in [26],



P σr21 (M1 (X̃)) ≥ (λ2 + r2 r3 )(1 − x) ≤ C exp Cr1 − c(λ2 + r2 r3 )x2 ∧ x ,

for any x > 0.

Similar results also hold for σmin (M2 (X̃)) and σmin (M3 (X̃)). Let x = 1/2, note that λ ≥
Cgap p1/2 , we have
σmin (Mk (X̃)) ≥

λ
21/2

,

k = 1, 2, 3

(56)

with probability at least 1 − C exp(−cp) for large enough constant Cgap > 0. Additionally,
by Lemma 5, we have
max

Vk ∈Opk ,rk
k=1,2,3

Z ×1 V1> ×2 V2> ×3 V3>

with probability at least 1 − C exp(−cp).
2

F

√
≤ C pr,

(57)

2. In the following Steps 2 and 3, we temporarily ignore the randomness of Z and the definition
of Û1 , Û2 , Û3 as the the estimators with high likelihood values. Instead we only assume
X, Z ∈ Rp1 ×p2 ×p3 , X = S ×1 U1 ×2 U2 ×3 U3 , and Ûk ∈ Opk ,rk for k = 1, 2, 3 satisfy (55), (56),
> Û k . Our goal in Steps 2 is to show
and (57). By Lemma 1 in [26], k sin Θ(Ûk , Uk )kF = kUk⊥
k F

under such setting, one must has
√
k sin Θ(Ûk , Uk )kF =

>
kUk⊥
Ûk kF

≤C

pr
,
λ

k = 1, 2, 3.

(58)

To simplify our notations, we first perform spectral transformation on (Ûk )[1:rk ,:] and (Ûk )[(rk +1):pk ,:] .
To be specific, let (Ûk )[1:rk ,:] = Ūk Σ̄k V̄k> be the singular value decomposition, and (Ûk )[(rk +1):pk ,:] V̄k =
Q̄k R̄k be the full QR decomposition (so that Q̄k is a square orthogonal matrix). Then we
transform


(Ûk )[1:rk ,:]


(Ûk )[(rk +1):pk ,:]

⇒

Z

⇒

X

⇒


Ū >
 k

 

(Ûk )[1:rk ,:]





Ūk> (Ûk )[1:rk ,:] V̄k





Σ̄k



 =  ,
·
 V̄k = 
R̄k
(
Û
)
V̄
Q̄>
(Ûk )[(rk +1):p,:]
k [(rk +1):pk ,:] k
k





Ū
Ū
 ×2  2
 ×3  3
,
Q̄1
Q̄2
Q̄3






Ū1
Ū2
Ū3
 ×2 
 ×3 
.
X ×1 
Q̄1
Q̄2
Q̄3

Q̄>
k

Ū1
Z ×1 

We can check that (55), (56), (57) still hold after this transformation.
Suppose diag(R1 ) =
 
Σ̄1
(a1 , . . . , ar ). Since Σ̄1 is diagonal and R̄1 is upper diagonal,   is orthogonal, we must have
R̄1
p
p
all off-diagonal entries of R̄k are zero, and Σ̄1 = diag( 1 − a21 , . . . , 1 − a2r ). For convenient
we also denote
(0)

ai

=

q
q
q
(1)
(0)
(1)
(0)
(1)
1 − a2i , ai = ai , bj = 1 − b2j , bj = bj , ck = 1 − c2k , ck = ck .

(59)

Therefore, without loss of generality we can assume there exist real numbers 0 ≤ ai , bj , ck ≤ 1

3

such that

(0)
a
 1






 (1)
Û1 = 
a1










..

.

..

.


(0)
b

 1








(0) 

ar1 



 , Û2 = b(1)

 1









(1) 

ar1 



0p1 −2r1 ,r1



..

.

..

.


(0)
c

 1








(0) 

br2 



 , Û3 = c(1)

 1









(1) 

br2 



0p2 −2r2 ,r2


..

.

..

.





(0) 
cr3 


.





(1) 
cr3 


0p3 −2r3 ,r3
(60)

where 0pk −2rk ,rk represents the zero matrix with dimension (pk − 2rk )-by-rk . By the form of
U1 , U2 , U3 in (53), we must have
r1
r2
r3
X
X
X
max k sin Θ(Ûk , Uk )kF = max (
a2i )1/2 , (
b2i )1/2 , (
c2i )1/2

(

k=1,2,3

i=1

i=1

)
.

(61)

i=1

 P 1 2 1/2 Pr2 2 1/2 Pr3 2 1/2
≤
ai ) , ( i=1 bi ) , ( i=1 ci )
In order to show (58) we only need to prove that max ( ri=1
√
C pr.
Next, we decompose the noise tensor Z to the following eight pieces,
Z(t1 t2 t3 ) = Z[(t1 r1 +1):(t1 r1 +r1 ),(t2 r2 +1):(t2 r2 +r2 ),(t3 r3 +1):(t3 r3 +r3 )] ,

t1 , t2 , t3 ∈ {0, 1}.

(62)

√
By (57), we know kZ(t1 t2 t3 ) kF ≤ C pr, t1 , t2 , t3 ∈ {0, 1}. Based on the form of Û1 , Û2 , Û3
(60), we have


Y ×1 Û1 ×2 Û2 ×3 Û3

(0) (0) (0)
=X̃ijk ai bj ck

+


ijk

(100) (1) (0) (0)
Zijk ai bj ck

(110) (1) (1) (0)

(010) (0) (1) (0)

(001) (0) (0) (1)

+ Zijk ai bj ck + Zijk ai bj ck

(101) (1) (0) (1)

(011) (0) (1) (1)

(111) (1) (1) (1)

+ Zijk ai bj ck + Zijk ai bj ck + Zijk ai bj ck + Zijk ai bj ck
X
(0) (0) (0)
(t t t ) (t ) (t ) (t )
=X̃ijk ai bj ck +
Zijk1 2 3 ai 1 bj 2 ck 3 .
t1 ,t2 ,t3 ∈{0,1}
t1 ,t2 ,t3 are not all 0

4

Therefore,
(55)

0 ≤

=

Y ×1 Û1 ×2 Û2 ×3 Û3

r1X
,r2 ,r3
i,j,k=1

≤

r1X
,r2 ,r3

2
F

− kY ×1 U1 ×2 U2 ×3 U3 k2F
2


X̃ijk a(0) b(0) c(0) +
i
j
k

2
X̃ijk




(0) (0) (0) 2
ai bj ck

X

(t t t3 ) (t1 ) (t2 ) (t3 ) 
ai bj ck 


Zijk1 2

t1 ,t2 ,t3 ∈{0,1}
t1 ,t2 ,t3 are not all 0



2 
− X̃ijk



−1

i,j,k=1

X

+2

+7

r1X
,r2 ,r3

X̃ijk

i,j,k=1
t1 ,t2 ,t3 ∈{0,1}
t1 ,t2 ,t3 are not all 0
r1X
,r2 ,r3 
X

max


(0) (0) (0) 2 (t1 t2 t3 ) (t1 ) (t2 ) (t3 )
ai bj ck
Zijk
ai bj ck


(t t t3 ) (t1 ) (t2 ) (t3 ) 2
ai bj ck

Zijk1 2

i,j,k=1
t1 ,t2 ,t3 ∈{0,1}
t1 ,t2 ,t3 are not all 0
r1X
,r2 ,r3

2
≤
X̃ijk
(1 − a2i )(1 − b2j )(1
i,j,k=1
(
r1X
,r2 ,r3

+ 63



− c2k ) − 1

X̃ijk

q

t1 ,t2 ,t3 ∈{0,1}
t1 ,t2 ,t3 are not all 0 i,j,k=1
)
r1X
,r2 ,r3 

(t1 t2 t3 ) (t1 ) (t2 ) (t3 ) 2
Zijk
ai bj ck
.
i,j,k=1



1 − a2i

q
q
(t t t ) (t ) (t ) (t )
1 − b2j 1 − c2k Zijk1 2 3 ai 1 bj 2 ck 3 ,

By the inequality above, one of the following inequalities must hold for some t1 , t2 , t3 ∈ {0, 1}
and t1 , t2 , t3 are not all 0:
0≤

r1X
,r2 ,r3



2
(1 − a2i )(1 − b2j )(1 − c2k ) − 1
X̃ijk
(63)

i,j,k=1

q
q
q
(t t t ) (t ) (t ) (t )
2
2
+ 63X̃ijk 1 − ai 1 − bj 1 − c2k Zijk1 2 3 ai 1 bj 2 ck 3 ,
0≤

r1X
,r2 ,r3

r1X
,r2 ,r3 



(t t t ) (t ) (t ) (t ) 2
2
X̃ijk
(1 − a2i )(1 − b2j )(1 − c2k ) − 1 + 63
Zijk1 2 3 ai 1 bj 2 ck 3
.

i,j,k=1

(64)

i,j,k=1

3. Next we discuss in two different situations: (63) or (64) hold.
(a) When (63) holds, we first assume t1 = 1, t2 = t3 = 0 as the other situations follow

5

similarly. Then
0≤

r1X
,r2 ,r3

2
X̃ijk



(1 −

a2i )



− 1 + 63

(100)

|X̃ijk ||Zijk |ai

(since 1 ≤ ai , bj , ck ≤ 1)

i,j,k=1

i,j,k=1

≤−

r1X
,r2 ,r3

r1X
,r2 ,r3


2
X̃ijk
a2i + 63 

r1X
,r2 ,r3

1/2 

r1X
,r2 ,r3

2
X̃ijk
a2i 



Thus,

r1X
,r2 ,r3

(100)

(Zijk )2 

,

(Cauchy-Schwarz inequality)

i,j,k=1

i,j,k=1

i,j,k=1

1/2

2
X̃ijk
a2i ≤ C

i,j,k=1

r1X
,r2 ,r3

(100)

(Zijk )2 .

(65)

i,j,k=1

Additionally, since
rX
2 ,r3

2
X̃ijk
=




M1 (X̃)

j,k=1

[i,:] 2

2
≥ σmin
(M1 (X̃))

and (56), we have
r1
X
i=1

λ
a2i

2

2

≤

r1X
,r2 ,r3

2
a2i
X̃ijk

(61)

0≤

2
[(1
X̃ijk

−

a2i )(1

−

(100)

Zijk

(57)

= 63kZ(100) k2F ≤ Cpr,

i,j,k=1

i,j,k=1

which means k sin Θ(Û1 , U1 )kF =
r1X
,r2 ,r3

≤ 63

r1X
,r2 ,r3

b2j )

qP
r1

2
i=1 ai

− 1] + 63

√
≤ C pr/λ. On the other hand, by (63),

r1X
,r2 ,r3

(100)
|X̃ijk ||Zijk |ai

q
1 − a2i (1 − b2j )

i,j,k=1

i,j,k=1

(by Algorithmic-geometric inequality)

r1X
,r2 ,r3 
r1X
,r2 ,r3
632 (100) 2
2
2
a2i (1 − b2j ) +
[−b2j − a2i + a2i b2j ] +
X̃ijk
≤
X̃ijk
(Zijk ) (1 − a2i )(1 − b2j )
4
i,j,k=1
r1X
,r2 ,r3

≤−

i,j,k=1

i,j,k=1

2 2
X̃ijk
bj +

632
4

r
X

(100)

(56)

(Zijk )2 ≤ −

r2
λ2 X
632 (100) 2
b2j +
kZ
kF .
2
4
j=1

i,j,k

(66)
Therefore,
k sin Θ(Û2 , U2 )kF

v
uX
(57) √
u r2 2
bj ≤ CkZ(100) kF /λ ≤ C pr/λ.
= t

(61)

j=1

√
By symmetry, one can also show that k sin Θ(Û3 , U3 )kF ≤ C pr/λ. In summary, we
must have

√
max k sin Θ(Ûk , Uk )kF ≤ C

k=1,2,3

for some constant C > 0 when (63) holds.
6

pr
λ

(b) When (64) holds for some t1 , t2 , t3 ∈ {0, 1},
0≤

≤

r1X
,r2 ,r3
i,j,k=1
r1X
,r2 ,r3

2
X̃ijk
[(1

−

a2i )(1

−

b2j )(1

−

c2k )

− 1] + 63

r1X
,r2 ,r3

(t t t3 ) 2

(Zijk1 2

(t1 ) 2

) (ai

(t )

(t )

) (bj 2 )2 (ck 3 )2

i,j,k=1
2
X̃ijk
[(1 − a2i ) − 1] + 63

i,j,k=1
r1X
,r2 ,r3

≤−

r1X
,r2 ,r3

(t t t3 ) 2

(Zijk1 2

)

i,j,k=1
2
+ 63kZ(t1 t2 t3 ) k2F ≤
a2i X̃ijk

i,j,k

r1 ,r2 ,r3
λ2 X 2
ai + Cpr,
2
i,j,k=1

(61)

which means k sin Θ(Û1 , U1 )kF =

qP
r1

2
i=1 ai

√
≤ C pr/λ. One can similarly prove the

parallel results for k sin Θ(Û2 , U2 )kF and k sin Θ(Û3 , U3 )kF . In summary, we also have
√
max k sin Θ(Ûk , Uk )kF ≤ C

k=1,2,3

pr
λ

for some constant C > 0 when (64) holds.
To sum up, we have the derived perturbation bound: under (55), (56), (57), one must have
√
k sin Θ(Ûk , Uk )kF ≤ C pr/λ.
4. Next we consider the recovery loss for X̂. Similarly as Steps 2-3, we temporarily ignore the
randomness of Z, and the definition of Û1 , Û2 , Û3 as the estimators with high likelihood values
in this step. We aim to prove
X̂ − X

F

√
≤ C pr,

under the assumptions of (55), (56), and (57). First, without loss of generality we can assume
Û1 , Û2 , Û3 have the simple form (60). Based on the structure of Û1 , Û2 , Û3 , we know

(0)
(1) (0)
(a )2
a1 a1
 1

..
..

.
.
0r1 ,p1 −2r1


(0)
(1) (0)

(ar1 )2
ar1 ar1

 (1) (0)
(1)
PÛ1 = 
(a1 )2
a1 a1

..
..

.
.
0r1 ,p1 −2r1



(1) 2
(1) (0)

(ar1 )
ar ar1

0p1 −2r1 ,r1
0p1 −2r1 ,r1
0p1 −2r1 ,p1 −2r1










,









while PÛ2 and PÛ3 can be written in similar forms. We have the following decomposition for

7

kX̂ − XkF ,
X̂ − X

F

= Y ×1 PÛ1 ×2 PÛ2 ×3 PÛ3 − S ×1 U1 ×2 U2 ×3 U3

≤ Z ×1 PÛ1 ×2 PÛ2 ×3 PÛ3

F

F

+ X ×1 PÛ1 ×2 PÛ2 ×3 PÛ3 − S ×1 U1 ×2 U2 ×3 U3

F

√

≤C pr + S ×1 (PÛ1 U1> ) ×2 (PÛ2 U2> ) ×3 (PÛ3 U3> ) − S ×1 U1 ×2 U2 ×3 U3
√

≤C pr + X̃

×1 (PÛ1 U1> )

×2 (PÛ2 U2> )

×3 (PÛ3 U3> )

− X̃ ×1 U1 ×2 U2 ×3 U3

F
F

+ Z[1:r,1:r,1:r] ×1 (PÛ1 U1> ) ×2 (PÛ2 U2> ) ×3 (PÛ3 U3> ) + Z[1:r,1:r,1:r] ×1 U1 ×2 U2 ×3 U3
√

≤C pr + X̃ ×1 (PÛ1 U1> ) ×2 (PÛ2 U2> ) ×3 (PÛ3 U3> ) − X̃ ×1 U1 ×2 U2 ×3 U3

F

F

.
(67)

Based on the form of U1 , U2 , U3 , Û1 , Û2 , Û3 , we have
2

X̃ ×1 (PÛ1 U1> ) ×2 (PÛ2 U2> ) ×3 (PÛ3 U3> ) − X̃ ×1 U1 ×2 U2 ×3 U3
F

2
r1X
,r2 ,r3
r1X
,r2 ,r3
2


X
(0) (0) (0)
(0) (0) (0) (t ) (t ) (t ) 2
2
2
ai bj ck
=
X̃ijk
−1 +
ai bj ck ai 1 bj 2 ck 0
X̃ijk
i,j,k=1
t1 ,t2 ,t3 ∈{0,1}
t1 ,t2 ,t3 are not all 0

i,j,k=1

=

r1X
,r2 ,r3





 
(0)
(1)
(0)
(1)
(0) (0) (0) 2
(0)
(1)
2
ai bj ck
(ai )2 + (ai )2 (bj )2 + (bj )2 (ck )2 + (ck )2
X̃ijk

i,j,k=1
r1X
,r2 ,r3

+

 


(0) (0) (0) 2
2
X̃ijk
−2 ai bj ck
+1 .

i,j,k=1
(0)

(1)

(0)

(1)

(0)

(1)

Recall the actual values of ai , ai , bj , bj , ck , ck in (59), we further have
X̃ ×1 (PÛ1 U1> ) ×2 (PÛ2 U2> ) ×3 (PÛ3 U3> ) − X̃ ×1 U1 ×2 U2 ×3 U3
 

r1X
,r2 ,r3
r1X
,r2 ,r3



(0) (0) (0) 2
(0) (0) (0) 2
2
2
X̃ijk ai bj ck
+
X̃ijk −2 ai bj ck
+1
=

i,j,k=1
r1X
,r2 ,r3

i,j,k=1

2
F

(68)


2
X̃ijk
1 − (1 − a2i )(1 − b2j )(1 − c2k ) .

i,j,k=1

By the analysis in Step 2, we know under (55), (56), (57), at least one of (63) and (64) must
hold for some (t1 , t2 , t3 ) ∈ {0, 1}3 \{(0, 0, 0)}. Again, we discuss in two different situations to
show no matter which of (63) or (63) happen, we must have
r1X
,r2 ,r3



2
X̃ijk
1 − (1 − a2i )(1 − b2j )(1 − c2k ) ≤ Cpr.

i,j,k=1

8

(69)

(a) When (63) holds, we again assume t1 = 1, t2 = t3 = 0 as the other situations follow
similarly. Particularly, we have shown in Step 2 (a), (65) and (66),
r1X
,r2 ,r3

2
X̃ijk
a2i

≤C

i,j,k=1

Clearly,

r1X
,r2 ,r3

r1X
,r2 ,r3

(100)
(Zijk )2 ,

i,j,k=1

Pr1 ,r2 ,r3
i,j,k=1

2 c2 ≤ C
X̃ijk
k
r1X
,r2 ,r3

=

≤

i,j,k=1
r1X
,r2 ,r3
i,j,k=1
r1X
,r2 ,r3

Pr1 ,r2 ,r3
i,j,k=1

2 2
X̃ijk
bj

≤C

r1X
,r2 ,r3

i,j,k=1

(100)

(Zijk )2

i,j,k=1

(100)

(Zijk )2 can be derived by symmetry. Then,


2
X̃ijk
1 − (1 − a2i )(1 − b2j )(1 − c2k )
2
X̃ijk
a2i + b2j + c2k − a2i b2j − a2i c2k − b2j c2k + a2i b2j c2k

2
X̃ijk
(a2i + b2j + c2k )

i,j,k=1
r1X
,r2 ,r3

≤C



(since 0 ≤ ai , bj , ck ≤ 1)

(100)

(Zijk )2 ≤ Cpr.

i,j,k=1

(b) When (64) holds, one has
r1X
,r2 ,r3

r1X
,r2 ,r3 


(t t t ) (t ) (t ) (t ) 2
2
Zijk1 2 3 ai 1 bj 2 ck 3
1 − (1 − a2i )(1 − b2j )(1 − c2k ) ≤ 63
X̃ijk

i,j,k=1
r1X
,r2 ,r3 

≤C

ijk=1
(t t t3 )

Zijk1 2

2

= CkZ(t1 t2 t3 ) kF ≤ Cpr.

i,j,k=1

In summary of Cases (a)(b), we must have (69). Combining (68), (67), and (69), we have
shown
√
kX̂ − XkF ≤ C pr,

under (55), (56), (57).

(70)

5. We finalize the proof for Theorem 2 in this step. We let Q = {(55), (56), (57) all hold}. By
√
Step 1, P (Q) ≥ 1−C exp(−cp); by Steps 2-4, one has k sin Θ(Ûk , Uk )kF ≤ C pr/λ, k = 1, 2, 3,
√
and kX̂ − XkF ≤ C pr under Q. The rest of the proof is essentially the same as the Step 4
in the proof of Theorem 1.
Since X̂ is a projection of Y by definition, so kX̂kF ≤ kYkF ≤ kXkF + kZkF . Then we have
the following upper bound for 4-th moment of recovery error,


EkX̂ − Xk4F ≤ C EkX̂k4F + kXk4F ≤ CkXk4F + CEkZk4F
 2
≤C exp(c0 p) + CE χ2p3 = C exp(c0 p) + Cp6 .
9

The we have the following upper bound for the risk of X̂,
EkX̂ −

Xk2F

= EkX̂ −

Xk2F 1Q

+ EkX̂ −

Xk2F 1Qc

q
= Cpr + EkX̂ − Xk4F E1Qc

≤Cpr + C exp ((c0 − c)p) + Cp6 exp(−cp).
Thus, one can select c0 < c to ensure that
EkX̂ − Xk2F ≤ Cpr ≤ p1 r1 + p2 r2 + p3 r3 .
Additionally, when σmin (Mk (X)) ≥ λ, we have kXk2F = kMk (X)k2F ≥ rk λ2 for k = 1, 2, 3,
which implies kXk2F ≥ Crλ2 . Thus we also have
E

kX̂ − Xk2F
p1 + p2 + p3
≤
.
2
λ2
kXkF

Now we consider the Frobenius sin θ norm risk for Ûk . Since sin Θ(Ûk , Uk ) is a rk -by-rk matrix
with spectral norm no more than 1, definition k sin Θ(Ûk , Uk )k2F ≤ rk ≤ r. Therefore, one has
Ek sin Θ(Ûk , Uk )kF = Ek sin Θ(Ûk , Uk )kF 1Q + Ek sin Θ(Ûk , Uk )kF 1Qc
√
√
q
pr
pr p
2
=C
+ Ek sin Θ(Ûk , Uk )kF · E1Qc ≤ C
+ r · C exp(−cp)
λ
λ
By the definition of λ, we know λ = σrk (Mk (X)) ≤

kXkF
√
rk

p)
√ 0 , so we can select c0 > 0
≤ C exp(c
rk

small enough to ensure that
√

p
p
pr
pr2
≥
≥ c r · C exp(−cp).
λ
C exp(c0 p)
√

pr
λ .

This means Ek sin Θ(Ûk , Uk )kF ≤ C
−1/q
Erk

sin Θ(Ûk , Uk )

q

≤

Finally, for any 1 ≤ q ≤ 2, we have

−1/2
Erk

sin Θ(Ûk , Uk )

To sum up, we have finished the proof for this theorem.

A.2

F

√
√
p
pk
≤C
≤C
.
λ
λ



Proof of Proposition 1

For convenient, we introduce the following notations: mk = |C ∩ Dk |, u(k) = (1C )Dk , U (k) =
u(k) /ku(k) k

2

is the normalized vector of

u(k) ,

(k)
U⊥

10

∈R

N
6

×



N2
−1
36



is the orthogonal complement

of U (k) , Ãk = Mk (2 · A[D1 ,D2 ,D3 ] − 1|D1 |×|D2 |×|D3 | ), for k = 1, 2, 3. Without loss of generality and
for convenience of the presentation, we



1,



(Ã1 )i,p3 (j−1)+k =
1,




 −1,

assume N is a multiple of 6. Based on the statement,


(1) (2) (3)
w.p 1, if ui , uj , uk
= (1, 1, 1);


(1) (2) (3)
(71)
w.p 1/2, if ui , uj , uk
6= (1, 1, 1);


(1) (2) (3)
w.p 1.2, if ui , uj , uk
6= (1, 1, 1).

A2 and A3 have the similar form. Therefore, Ã1 are all 1 in the block of (D1 ∩ C) × ((D2 ∩
C) ⊗ (D3 ∩ C)), and are with i.i.d. Rademacher entries outside the block. Since C is uniformly
randomly selected from V1 , |V1 | = N/6, |Dk | = N/6, |C| = κN , we know m1 = |D1 ∩ C|, m2 =
|D1 ∩ C|, m3 = |D1 ∩ C| satisfy hypergeometric distribution with parameter (κN , N/2, N/6).
Based on the concentration inequality of hypergeometric distribution (Theorem 1 in [64]),
κN
κN
≤ mk = |Dk ∩ C| ≤
,
4
2

k = 1, 2, 3

(72)

with probability at least 1 − C exp(−cκN ). Now the rest of the proof is similar to Theorem 3 in
[26]. By (71), we have




(U (1) )> Ã1 ∈ RN

2 /36

,




  = √m1 , j ∈ ((D2 ∩ C) ⊗ (D3 ∩ C));
(U (1) )> A1
j
∼ √W
otherwise,
m1 ,

where W has the same distribution as the sum of m1 i.i.d. Rademacher random variables.
Conditioning on C satisfying (72), similarly as the derivation for Equation (1.15) in the Appendix
of [26], we can derive

 N2 m m m
1 2 3
+
σ12 (U (1) )> A1 ≥
36
4

with probability at least 1 − C exp(−cN );

N 2 m1 m2 m3
+
with probability at least 1 − C exp(−cN );
36
8
√
(1)
k(u⊥ )> A1 P(U (1) )> A1 k ≤ C N with probability at least 1 − C exp(−cN ).
σ22 (A1 ) ≤

(73)
(74)
(75)

Under the circumstance that (72), (73), (73), and (73) all hold, by Proposition 1 in [26], we
have
(1)

σ2 (A1 )k(u⊥ )> A1 P(U (1) )> A1 k

≤
σ12 (U (1) )> A1 − σ22 (A1 )
p
3/2
(73)(73)(73)
N (N 2 + m1 m2 m3 ) (72) N 3/2 + N 1/2 κN
≤ C
.
≤
C
m1 m2 m3
κ3N

(1)
k sin Θ(û>
)k
1 ,U

11

√
Note that lim inf N →∞ κN / N = ∞, limN →∞ P ((72), (73), (73), and (73) all hold) = 1, we
have
d

>
(1)
)k → 0,
sin Θ(û>
1 , (1C )D1 ) = k sin Θ(û1 , U

as N → ∞.

The proofs for k = 2, 3 essentially follow. Therefore, we have finished the proof of this proposition.

B



Appendix: Technical Lemmas

We collect all technical lemmas that has been used in the theoretical proofs throughout the
paper in this section.
The following lemma shows the equivalence between two widely considered Schatten q-norm
distances for singular subspaces.
Lemma 3. For any U1 , U2 ∈ Op,r and all 1 ≤ q ≤ +∞,
1
kU1 U1> − U2 U2> kq ≤ k sin Θ(U1 , U2 )kq ≤ kU1 U1> − U2 U2> kq .
4
We will use the following properties of tensor algebra in the technical analysis of this paper.
Lemma 4 (Properties in Tensor Algebra).
• Suppose X ∈ Rp1 ×p2 ×p3 , Uk ∈ Rpk ×rk for k = 1, 2, 3. Then we have the following identity
related to tensor matricizations,


>
>
Mk X ×k+1 Uk+1
×k+2 Uk+2
= Mk (X)(Uk+1 ⊗ Uk+2 ),

k = 1, 2, 3.

(76)

• Suppose we further have Ũk ∈ Rpk ×r̃k for k = 1, 2, 3, then


Ũ1 ⊗ Ũ2

>

= (Ũ1> ) ⊗ (Ũ2> ),



Ũ2 ⊗ Ũ3

kU2 ⊗ U3 k = kU2 k · kU3 k,

>


 

(U2 ⊗ U3 ) = Ũ2> U2 ⊗ Ũ3> U3 .

kU2 ⊗ U3 kF = kU2 kF · kU3 kF ,

σmin (U2 ⊗ U3 ) = σmin (U2 )σmin (U3 ).

12

(77)

(78)

• (Properties related to projections) Suppose U2 ∈ Op2 ,r2 , U3 ∈ Op3 ,r3 , and U2⊥ ∈ Op2 ,p2 −r2 , U3⊥ ∈
Op3 ,p3 −r3 are their orthogonal complement, respectively. Then PU2 ⊗U3 = PU2 ⊗ PU3 , and
we have the following decomposition
Ip2 p3 =PIp2 ⊗U3 + PIp3 ⊗U3⊥ = PU2 ⊗Ip3 + PU3⊥ ⊗Ip2

(79)

=PU2 ⊗U3 + PU2⊥ ⊗U3 + PU2 ⊗U3⊥ + PU2⊥ ⊗U3⊥ .
The following lemma characterizes the maximum of norms for i.i.d. Gaussian tensors after
any projections.
iid

Lemma 5. For i.i.d. Gaussian tensor Z ∈ Rp1 ×p2 ×p3 , Z ∼ N (0, 1), we have the following tail
bound for the projections,


P

max

V2 ∈Rp2 ×r2 ,V3 ∈Rp3 ×r3
kV2 k≤1,kV3 k≤1




M1 Z ×2 V2> ×3 V3>

√
√
√
√
√

≥ C p1 + C r2 r3 + C 1 + t ( p2 r2 + p3 r3 )



≤C exp(−Ct(p2 r2 + p3 r3 ))


for any t > 0. Similar results also hold for M2 Z ×1 V1> ×3 V3> and M3 Z ×1 V1> ×2 V2> .
Meanwhile, there exists uniform C > 0 such that


P

max

∈Rp×r

V1 ,V2 ,V3
max{kV1 k,kV2 k,kV3 k}≤1



Z ×1 V1> ×2 V2> ×3 V3>

2
F


≥ Cr1 r2 r3 + C(1 + t)(p1 r1 + p2 r2 + p3 r3 )

≤ exp (−Ct(p1 r1 + p2 r2 + p3 r3 ))
(80)
for any t > 0.
In the perturbation bound analysis in this paper, we also need the following technical result
to bound the spectral and Frobenius norm for the projections.
Lemma 6. Suppose X, Z ∈ Rp1 ×p2 , rank(X) = r. If the singular value decomposition of X and
Y are written as
h

Y = X + Z = Û Σ̂V̂ > = Û1 Û2

13

i


Σ̂1
·


Σ̂2

h
 · V̂ >
1

V̂2>

i

,

where Û1 ∈ Op1 ,r , V̂1 ∈ Op2 ,r correspond to the leading r left and right singular vectors; and
Û2 ∈ Op1 ,p2 −r , V̂1 ∈ Op2 ,p2 −r correspond to their orthonormal complement. Then
PÛ2 X ≤ 2kZk,

PÛ2 X

F

 √
≤ min 2 rkZk, 2kZkF .

The following lemma provides a detailed analysis for ε-net for the class of regular matrices
under various norms and for the low-rank matrices under spectral norm.
Lemma 7 (ε-net for Regular and Low-rank Matrices).
• Suppose k · k• is any matrix norm, Xp1 ,p2 = {X ∈ Rp1 ×p2 : kXk ≤ 1} is the unit ball
around the center in k · k• norm. Then there exists an ε-net X̄p1 ,p2 in k · k• norm with
cardinality at most ((2 + ε)/ε)p1 p2 for Xp1 ,p2 . To be specific, there exists X (1) , . . . , X (N )
with N ≤ ((2+ε)/ε)p1 p2 , such that for all X ∈ Xp1 ,p2 , there exists i ∈ {1, . . . , N } satisfying
kX (i) − Xk ≤ ε.
• Let Xp1 ,p2 ,r = {X ∈ Rp1 ×p2 : rank(X) ≤ r, kXk ≤ 1} be the class of low-rank matrices
under spectral norm. Then there exists an ε-net X̄r for Xp1 ,p2 ,r with cardinality at most
((4 + ε)/ε)(p1 +p2 )r . Specifically, there exists X (1) , . . . , X (N ) with N ≤ ((4 + ε)/ε)(p1 +p2 )r ,
such that for all X ∈ Xp1 ,p2 ,r , there exists i ∈ {1, . . . , N } satisfying kX (i) − Xk ≤ ε.
The next lemma characterizes the tail probability for i.i.d. Gaussian vector after multiplication of any fixed matrix.
iid

Lemma 8. Suppose u ∈ Rp such that x ∼ N (0, 1), A ∈ Rp×r is a fixed matrix, then

√
P kAuk22 − kAk2F ≤ −2kA> AkF t ≤ exp(−x);


√
P kAuk22 − kAk2F ≥ 2kA> AkF t + 2kAk2 t ≤ exp(−x).

C
C.1

Proof of Technical Lemmas
Proof of Lemma 1

Without loss of generality, assume that p ≡ 0(mod 2). Hereafter, set N = 3p and κN = 20k with
k = bp(1−τ )/2 c. Our main technique is based on a reduction scheme which maps any adjacency
14

tensor A ∈ {0, 1}N ×N ×N to a random tensor Y ∈ Rp×p×p in O(N 3 ) number of flops. The
technique was invented in [49], adapted from a bottom-left trick in [45]. Some other related
methods can be found in [50] and [47]. For the completeness and readability of our paper, we
provide a detailed application of this technique to the tensor settings.
To this end, for any M ≥ 3 and 0 < µ ≤
ξ + := (Z + µ)1(|Z| ≤ M )

1
2M ,

define two random variables

and ξ − := (Z̃ − µ)1(|Z̃| ≤ M )

where Z and Z̃ denote independent standard normal random variables. The randomized mapping from A ∈ {0, 1}N ×N ×N to a random matrix Y ∈ Rp×p×p is essentially one step of Gaussianization. For simplicity, denote V1 := {1, 2, . . . , p2 } ∪ { 3p
2 + 1, . . . , 2p},
V2 :=

p

5p
+ 1, . . . , p ∪ 2p + 1, . . . ,
,
2
2

and
 5p

3p
∪
+ 1, . . . , 3p .
V3 := p + 1, . . . ,
2
2
Therefore, V1 , V2 , V3 are disjoint and V1 ∪ V2 ∪ V3 = [N ]. Given an adjacency tensor A ∈
{0, 1}N ×N ×N , let A0 = AV1 ,V2 ,V3 ∈ Rp×p×p be a corner block of A. Conditioned on A0 , we
generate a random tensor Y ∈ Rp×p×p such that

+
Ya,b,c = 1 − (A0 )a,b,c Ξ−
a,b,c + (A0 )a,b,c Ξa,b,c ,

∀a, b, c ∈ [p]

where Ξ− ∈ Rp×p×p has i.i.d. entries with the same distribution as ξ − and Ξ+ ∈ Rp×p×p has
i.i.d. entries with the same distribution as ξ + . Clearly, this process defines a deterministic map
for any fixed Ξ− , Ξ+ ∈ Rp×p×p
T : {0, 1}N ×N ×N × Rp×p×p × Rp×p×p 7→ Rp×p×p
(A, Ξ− , Ξ+ ) 7→ Y.
Let L(X) denote the law of a random tensor X. The total variation distance between two
probability distributions P1 and P2 is denoted by dTV (P1 , P2 ). The following lemma is analogous
to [49, Lemma 2] and the proof is skipped here.
Lemma 9. Let M ≥ 4, µ ≤

1
2M ,

and η be a Bernoulli random variable. Suppose ξ is a random

variable such that (ξ|η = 1) = ξ + and (ξ|η = 0) = ξ − .
15


2
(1) If P(η = 1) = 1, then dTV L(ξ), N (µ, 1) ≤ e(1−M )/2 ;
(2) If P(η = 0) = P(η = 1) = 21 , then dTV (L(ξ), N (0, 1)) ≤ e−M
Our next step is to show that (by choosing M =

√

2 /2

.

8 log 3p and µ = (2M )−1 ), the law of

Y = T (A, Ξ− , Ξ+ ) is asymptotically equivalent to a mixture over {PX : X ∈ M0 (p, k, r, λ)}
for λ =

3(1−τ )/4
p√
2 8 log 3p

if G ∼ H0 . On the other hand, if G ∼ H1 , the law of Y = T (A, Ξ− , Ξ+ ) is

asymptotically equivalent to a mixture over {PX : X ∈ M1 (p, k, r, λ)}. For an adjacency tensor
A ∈ {0, 1}N ×N ×N , we have T (A, Ξ− , Ξ+ ) ∈ Rp×p×p . Recall that Y ∈ Rp1 ×p2 ×p3 and we define
an embedding ` : Rp×p×p 7→ Rp1 ×p2 ×p3 ,


Aijk
`(A)ijk =

0

if (i, j, k) ∈ [p] × [p] × [p];
(81)
otherwise.

Lemma 10 is similar to [49, Lemma 4]. We postpone the proof of Lemma 10 to the Appendix.
Lemma 10. Let A ∈ RN ×N ×N be the adjacency tensor of a hypergraph G sampled from either
√
1
H0 or H1 and Y = ` ◦ T (A, Ξ− , Ξ+ ). Suppose that M = 8 log N and µ = 2M
. For each
i = 0, 1, if G ∼ Hi , there exists a prior distribution πi on Mi (p, k, r, λ) with λ =
that

3/4(1−τ )
p√
2 8 log 3p

such

√



dTV L(Y), Pπi ≤
where Pπi =

R

2.5k
e
+ 6k 0.86
.
27N

Mi (p,k,r,λ) PX (·)πi (dX).

Now, on the contradictory, suppose that the claim of Lemma 1 does not hold. It means that
there exists a sequence of polynomial-time tests {φpt } with a sub-sequence (pt )∞
t=1 of positive
integers such that
lim Rp,r,λ (φpt ) = lim

t→∞

t→∞

n

sup


PX φpt (Y) = 1 +

X∈M0 (p,k,r,λ)

o 1

PX φpt (Y) = 0
< .
2
X∈M1 (p,k,r,λ)
sup


Define the test ψNt (A) = φpt ` ◦ T (A, Ξ− , Ξ+ ) and we obtain a sequence of polynomial-time
tests {ψNt } for problem (16) with Nt = 3pt for t = 1, . . . , ∞. It suffices to compute


RNt ,κNt (ψNt ) = PH0 ψNt (A) = 1 + PH1 ψNt (A) = 0
(1−τ )/2

with κNt = 20bpt

c. Note that limt→∞

log κNt
√
log Nt

under H0 ,


PH0 ψNt (A) = 1 − Pπ0 φpt (Y) = 1
16

≤ 1−τ . By definition of dTV and Lemma 10,




= PH0 φpt ` ◦ T (A, Ξ− , Ξ+ ) = 1 − Pπ0 φpt (Y) = 1
√


e
−
+
≤dTV L ` ◦ T (A, Ξ , Ξ ) , Pπ0 ≤
+ 6kt (0.86)2.5kt
27Nt
(1−τ )/2

where kt = bpt

c and we used the fact that the mixture Pπ0 over M0 (p, k, r, λ) is also a

mixture over Fp,r (λ). In a similar fashion,
√



PH1 ψNt (A) = 0 − Pπ1 φpt (Y) = 0

≤

e
+ 6kt (0.86)2.5kt
27Nt

As a result,


RNt ,κNt (ψNt ) = PH0 ψNt (A) = 1 + PH1 ψNt (A) = 0
√


2 e
≤ Pπ0 φpt (Y) = 1 + Pπ1 φpt (Y) = 0 +
+ 12kt (0.86)2.5kt
27Nt
√


2 e
+ 12kt (0.86)2.5kt
≤
sup
PX φpt (Y) = 1 +
sup
PX φpt (Y) = 0 +
27N
t
X∈M0 (p,k,r,λ)
X∈M1 (p,k,r,λ)
Therefore,
1
lim RNt ,κNt (ψNt ) < ,
t→∞
2
which contradicts the hypothesis H(τ ).

C.2



Proof of Lemma 3

Let σ1 ≥ σ2 ≥ . . . ≥ σr denote the singular values of U1> U2 . It is easy to check that the singular
p
p
p
> U are
values of U1⊥
1 − σr2 ≥ . . . ≥ 1 − σ22 ≥ 1 − σ12 , in view of the fact
2
>
>
U2 )> (U1⊥
U2 ) = U2> U2 = Ir .
(U1> U2 )> (U1> U2 ) + (U1⊥

Recall that for all 1 ≤ q ≤ +∞,
k sin Θ(U1 , U2 )kq =

r
X

sin(cos−1 σi )

q 1/q

=

r
X

1 − σi2

q/2 1/q

.

i=1

i=1

The following fact is straightforward:
kU2 U2>

−

U1 U1> kq

≥

>
kU1⊥
U2 U2> kq

=

>
kU1⊥
U2 kq

=

r
X

1 − σi2

q/2 1/q

i=1

which concludes kU2 U2> − U1 U1> kq ≥ k sin Θ(U1 , U2 )kq . On the other hand,
kU2 U2> − U1 U1> kq ≤kPU1 (U2 U2> − U1 U1> )PU1 kq + kPU1 (U2 U2> )PU⊥1 kq
17

+kPU⊥1 (U2 U2> )PU1 kq + kPU⊥1 (U2 U2> )PU⊥1 kq

>
>
≤ U1 U1> U2 U2> U1 − Ir U1> q + kU2> U1⊥ kq + kU1⊥
U2 kq + kU1⊥
U2 U2> U1⊥ kq
r
X
q/2 1/q
≤ 4k sin Θ(U1 , U2 )kq
≤4
1 − σi2
i=1

where we used the fact 1 − σi2 ≤

C.3

q
1 − σi2 for all 1 ≤ i ≤ r.



Proof of Lemma 4


> ×
>
• First, we shall note that both Mk X ×k+1 Uk+1
k+2 Uk+2 and Mk (X)(Uk+1 ⊗ Uk+2 )
are of dimension pk -by-(rk+1 rk+2 ). To prove they are equal, we just need to compare
each of their entries. We focus on k = 1 as the k = 2, 3 essentially follows. For any
1 ≤ i1 ≤ p1 , 1 ≤ i2 ≤ r2 , 1 ≤ i3 ≤ r3 , one has
h

i
M1 X ×2 U2> ×3 U3>

i1 ,(i2 −1)r3 +i3

=

=

p2
X

p3
X

j2 =1 j3 =1
p2 X
p3
X



= X ×2 U2> ×3 U3>

i1 ,i2 ,i3

Xi1 ,j2 ,j3 (U2 )j2 ,i2 (U3 )j3 ,i3
(M1 (X))i1 ,(j2 −1)p3 +j3 · (U2 ⊗ U3 )(j2 −1)p3 +j3 ,(i2 −1)r3 +i3

j2 =1 j3 =1

= (M1 (X) · (U2 ⊗ U3 ))i1 ,(i2 −1)r3 +i3 .
This shows (76).
• The proof for (77) is essentially the same as (76) as we only need to check each entries of
the terms in (77) are equal. For (78), let
U2 =

X

>
σ2i · α2i β2i
,

i

U3 =

X

>
σ3j · α3j β3j

j

be the singular value decompositions. Then it is not hard to see the singular value decomposition of U2 ⊗ U3 can be written as
U2 ⊗ U3 =

X

σi σj · (αi2 ⊗ αj3 )(βi2 ⊗ βj3 )> ,

i,j

so that the singular values of U⊗ U3 are {σi · σj }. Then

 

kU2 ⊗ U3 k = max σi σj = max σi · max σj = kU2 k · kU3 k,
i,j

i

18

j

kU2 ⊗ U3 k2F =

X

σi2 σj2 =

X

i,j


! 
X
σi2 · 
σj2  = kU2 k2F · kU3 k2F ,

i

j

and

 

σmin (U2 ⊗ U3 ) = min σi σj = min σi · min σj = σmin (U2 ) · σmin (U3 ).
i,j

j

i

•
PIp2 ⊗U3 + PIp2 ⊗U3⊥ = (Ip2 ⊗ U3 )(Ip2 ⊗ U3 )> + (Ip2 ⊗ U3⊥ )(Ip2 ⊗ U3⊥ )>
(77)

>
>
)
) = Ip2 ⊗ (U3 U3> + U3⊥ U3⊥
= (Ip2 Ip>2 ) ⊗ (U3 U3> ) + (Ip2 Ip>2 ) ⊗ (U3⊥ U3⊥

=Ip2 ⊗ Ip3 = Ip2 p3 .
The other identity can be shown similarly.

C.4



Proof of Lemma 5

The key idea for the proof of this lemma is via ε-net. By Lemma 7, for k = 1, 2, 3, there exist
(1)

(Nk )

ε-nets: Vk , . . . , Vk

for {Vk ∈ Rpk ×rk : kVk k ≤ 1}, |Nk | ≤ ((4 + ε)/ε)pk rk , such that
(j)

For any V ∈ Rpk ×rk satisfying kV k ≤ 1, there exists Vk
(i)

For fixed V2

(j)

such that kVk

− V k ≤ ε.

(j)

and V3 , we consider
(ij)

Z1



(i)
(j)
= M1 Z ×2 (V2 )> ×3 (V3 )> ∈ Rp1 ×(r2 r3 ) .

 
 

(ij)
(i)> (i)
(j)> (j)
Clearly, each row of Z1 follows a joint Gaussian distribution: N 0, V2 V2
⊗ V3 V3
,

 

(i)> (i)
(j)> (j)
and V2 V2
⊗ V3 V3
≤ 1. Then by random matrix theory (e.g. [62]),


√
√
(ij)
P kZ1 k ≤ p1 + r2 r3 + t ≥ 1 − 2 exp(−t2 /2).
Then we further have


√
√
(ij)
P max kZ1 k ≤ p1 + r2 r3 + x ≥ 1 − 2((4 + ε)/ε)p2 r2 +p3 r3 exp(−x2 /2),
i,j

for all x > 0. Now, we assume
V2∗ , V3∗ =

M=

arg max
V2 ∈Rp2 ×r2 ,V3 ∈Rp3 ×r3
kV2 k≤1,kV3 k≤1

max

V2 ∈Rp2 ×r2 ,V3 ∈Rp3 ×r3
kV2 k≤1,kV3 k≤1

19



M1 Z ×2 V2> ×3 V3> ,


M1 Z ×2 V2> ×3 V3> .

(82)

(i)

By definition of the ε-net, we can find 1 ≤ i ≤ N2 and 1 ≤ j ≤ N3 such that kV2 − V2∗ k ≤ ε
(i)

and kV3 − V3∗ k ≤ ε. In this case under (82),


M = M1 Z ×2 (V2∗ )> ×3 (V3∗ )>




(i)
(j)
(i)
(j)
≤ M1 Z ×2 (V2 )> ×3 (V3 )> + M1 Z ×2 (V ∗ − V2 )> ×3 (V3 )>


(j)
+ M1 Z ×2 (V2∗ )> ×3 (V3∗ − V3 )>
√
√
≤ p1 + r2 r3 + x + εM + εM,
Therefore, we have


√
√
p1 + r2 r3 + x
≥ 1 − 2((4 + ε)/ε)p2 r2 +p3 r3 exp(−x2 /2).
P M≤
1 − 2ε
By setting ε = 1/3, x2 = 2 log(13)(p2 r2 + p3 r3 )(1 + t) for some large constant C > 0, we have
proved the first part of the lemma.
The proof for the second part is similar. For any given Vk ∈ Rpk ×rk satisfying kV1 k, kV2 k, kV3 k ≤
1, we have kV1 ⊗ V2 ⊗ V3 k ≤ 1. By Lemma 8, we know
P



2

Z ×1 V1> ×2 V2> ×3 V3> − kV1 ⊗ V2 ⊗ V3 k2F
F
q

≥ 2 tk(V1> V1 ) ⊗ (V2> V2 ) ⊗ (V3> V3 ) + 2tkV1 ⊗ V2 ⊗ V3 k2 ≤ exp(−t).

Since kV1 ⊗ V2 ⊗ V3 k ≤ 1, kV1 ⊗ V2 ⊗ V3 k2F = kV1 k2F kV2 k2F kV3 k2F ≤ r1 r2 r3 , then
k(V1> V1 ) ⊗ (V2> V2 ) ⊗ (V3> V3 )k2F = kV1> V1 k2F kV2> V2 k2F kV3> V3 k2F
! r
! r
!
r1
2
3
X
X
X
4
4
4
=
σi (V1 )
σi (V1 )
σi (V1 ) ≤ r1 r2 r3 ,
i=1

i=1

i=1

we have for any fixed V1 , V2 , V3 and x > 0 that


2
√
>
>
>
≥ r1 r2 r3 + 2 r1 r2 r3 x + 2x ≤ exp(−x).
P
Z ×1 V1 ×2 V2 ×3 V3
F

√
By geometric inequality, 2 r1 r2 r3 x ≤ r1 r2 r3 + x, then we further have


2
>
>
>
P
Z ×1 V1 ×2 V2 ×3 V3
≥ 2r1 r2 r3 + 3x ≤ exp(−x).
F

The rest proof for this lemma is similar to the first part. By Lemma 7, one can find three ε-nets:
(1)

(Nk )

Vk , . . . , Vk

for {Vk ∈ Rpk ×rk : kVk k ≤ 1} such that |Nk | ≤ ((4 + 2ε)/ε)pk rk , k = 1, 2, 3. Then

20

by probability union bound,

P

max
(a)

V1

Z

(c)

(b)

×1 V1>

×2 V2>

×3 V3>

,V2 ,V3

2
F


≥ 2r1 r2 r3 + 3x
(83)

≤ exp(−x) · ((4 + ε)/ε)p1 r1 +p2 r2 +p3 r3 .
When the inequality above holds, we suppose
(V1∗ , V2∗ , V3∗ ) = arg max

Vk ∈Rpk ×rk

Z ×1 V1> ×2 V2> ×3 V3> ,

and T = Z ×1 (V1∗ )> ×2 (V2∗ )> ×3 (V3∗ )> .

kVk k≤1
(a)

(b)

(c)

Then we can find V1 , V2 , V3

in the corresponding ε-nets such that

(a)

(b)

kV1∗ − V1 k ≤ ε,

kV2∗ − V2 k ≤ ε,

(c)

kV3∗ − V3 k ≤ ε.

Then
T = Z ×1 (V1∗ )> ×2 (V2∗ )> ×3 (V3∗ )>
(a)

(b)

(c)

(a)

≤ Z ×1 (V1 )> ×2 (V2 )> ×3 (V3 )> + Z ×1 (V1
(a)

(b)

− V1∗ )> ×2 (V2∗ )> ×3 (V3∗ )>
(a)

(b)

(c)

+ Z ×1 (V1 )> ×2 (V2 − V2∗ )> ×3 V3> + Z ×1 (V1 )> ×2 (V2 )> ×3 (V3


(a)
(b)
(c)
≤2r1 r2 r3 + 3t + kV1∗ − V1 k + kV2∗ − V2 k + kV3∗ − V3 k · T,

− V3∗ )>

which implies T ≤ (2r1 r2 r3 + 3x)/(1 − 3ε) provided that ε < 1/3 and (83) holds. Let ε = 1/9,
x = (1 + t) log(37) · (p1 r1 + p2 r2 + p3 r3 ) for some large constant C > 0, by (83) again we have
P (T ≥ Cr1 r2 r3 + C(1 + t)(p1 r1 + p2 r2 + p3 r3 ))

(84)

≥ exp(−Ct(p1 r1 + p2 r2 + p3 r3 ))
for some uniform constant C > 0. thus we have finished the proof for (80).

C.5



Proof of Lemma 6
PÛ2 X ≤ PÛ2 (X + Z) + kZk = σr+1 (Y ) + kZk =
≤kY − Xk + kZk = 2kZk.

Since rank (PU2 X) ≤ rank(X) ≤ r, it is clear that
√
kPU2 XkF ≤ 2 rkZk;
21

min

X̃∈Rp1 ×p2
rank(X̃)≤r

kY − X̃k + kZk

meanwhile,
PÛ2 X

F

≤ PÛ2 (X + Z)
≤

min

X̃∈Rp1 ×p2
rank(X̃)≤r

which has proved this lemma.

C.6

F

+ kZkF =

pX
1 ∧p2

!1/2
σi2 (Y

)

+ kZkF

i=r+1

kY − X̃kF + kZkF ≤ kY − XkF + kZkF ≤ 2kZkF ,



Proof of Lemma 7

• We first consider the ε-net for Xp1 ,p2 . Note that Xp1 ,p2 is a convex set in Rp1 ×p2 , we
sequentially pick matrices from Xp1 ,p2 , say X (1) , X (2) , . . . satisfying the following criterion:
for each time t, the picked matrix satisfies mint0 ≤t kX (t) − X (t−1) k• ≥ ε, i.e., the distances
from X (t) to all the other selected matrices are at least ε. We stop the selection process
until it is not possible to select the next matrix satisfying such criterion.
Suppose now X (1) , . . . , X (N ) are all we have selected. Since it is not possible to select
another matrix from Xp1 ,p2 which meets the criterion, all matrices in Xp1 ,p2 must be within
ε of some selected matrix in {X (1) , . . . , X (N ) }, thus
(i)
Xp1 ,p2 ⊆ ∪N
i=1 B(X , ε).

Here B(X (i) , ε) = {X ∈ Rp1 ×p2 : kX − X (i) k• ≤ ε} is the closed ball with center X (i) and
radius ε, Therefore, {X (1) , . . . , X (N ) } is a ε-net.
On the other hand, for any 1 ≤ i < j ≤ N , kX (i) − X (j) k• ≥ ε, so
(i)
{X ∈ Rp1 ×p2 : kXk• ≤ 1 + ε/2} ⊇ ∪N
i=1 B(X , ε/2),

and B(X (i) , ε/2) ∩ B(X (j) , ε/2) contains at most one matrix for any 1 ≤ i < j ≤ N .
Therefore,
(1 + ε/2)p1 p2 vol(Xp1 ,p2 ) = vol({X ∈ Rp1 ×p2 : kXk• ≤ 1 + ε/2})
≤

N
X

(85)
(i)

p1 p2

vol(B , ε/2) = N (ε/2)

i=1

which implies N ≤ ((2 + ε)/ε)p1 p2 .
22

vol(Xp1 ,p2 ),

• By the first part of this lemma, there exist (ε/2)-nets X̄p1 ,r and X̄r,p2 for {kX ∈ Rp1 ×r :
kXk ≤ 1} and {kX ∈ Rr×p2 : kXk ≤ 1}, such that




4 + ε p1 r
4 + ε p2 r
X̄p1 ,r ≤
, X̄r,p2 ≤
.
ε
ε
Next, we argue that

F̄p1 ,p2 ,r := X · Y : X ∈ X̄p1 ,r , Y ∈ X̄r,p2
is an ε-net for Fp1 ,p2 ,r in the spectral norm. Actually for any X ∈ Fp1 ,p2 ,r , we can find
A, B such that X = A · B, A ∈ Rp1 ×r , kAk ≤ 1; B ∈ Rr×p2 , kBk ≤ 1. Then we can find
A∗ ∈ X̄p1 ,r and B ∗ ∈ X̄r,p2 such that kA−A∗ k ≤ ε/2, kB −B ∗ k ≤ ε/2, thus A∗ B ∗ ∈ F̄p1 ,p2 ,r
satisfies
kX − A∗ B ∗ k = k(AB − AB ∗ ) + (AB ∗ − A∗ B ∗ )k
≤kAk · kB − B ∗ k + kA − A∗ k · kB ∗ k ≤ 1 · ε/2 + 1 · ε/2 = ε.
Note that X̄p1 ,p2 ,r ≤ X̄p1 ,p2 ,r · X̄p1 ,p2 ,r ≤ ((4 + ε)/ε)r(p1 +p2 ) , this has finished the proof
of this lemma.

C.7



Proof of Lemma 8

Suppose A = U ΣV > is the singular value decomposition of A. Since U, V are orthogonal and
P
iid
2 2
u ∼ N (0, 1), kAuk22 has the same distribution as p∧n
i=1 σi (A) ui . By the exponential probability
for general chi-square distribution (Lemma 1 in [65]), we have
v


u p∧n
p∧n
p∧n
X
X
u X
P
σi2 (A)u2i −
σi2 (A) ≤ −2tt
σi4 (A) ≤ exp(−x);
i=1

i=1

i=1

v


u p∧n
p∧n
p∧n
X
X
u X
σi2 (A)u2i −
σi2 (A) ≥ 2tt
σi4 (A) + 2t max σi2 (A) ≤ exp(−x),
P
i=1

i=1

i=1

which has finished the proof for Lemma 5 since kA> Ak2F =

Pp∧n
i=1

σi4 (A), and kAk = max σi (A).



C.8

Proof of Lemma 10

Clearly, it suffices to prove the claim for i = 0, i.e., under H0 . Let G = (V, E) ∼ H0 with
V = {1, 2, . . . , N } and A denote its adjacency tensor, meaning that there is a clique of size
23

κN planted in the subset {1, 2, . . . , bN/2c}. Recall that N = 3p for an even integer p. The

vertices set of the planted clique is denoted by C ⊂ 1, . . . , 3p
with |C| = κN = 20k where
2
k = bp(1−τ )/2 c. Recall V1 , V2 , V3 and define
Cj := C ∩ Vj ,

j = 1, 2, 3

which represents the subsets of clique vertices in V1 , V2 , V3 . If Y = T (A, Ξ− , Ξ+ ) ∈ Rp×p×p , it is
clear that, under H0 , X = E(Y|C) is a sparse tensor with supports S1 (X) = C1 ⊂ [p/2], S2 (X) =
C2 −

p
2

⊂ [p/2] and S3 (X) = C3 − p ⊂ [p/2]. We show that the sizes of S1 (X), S2 (X), S3 (X) are

lower bounded by k with high probability.
Lemma 11. There exists an event E on which min{|S1 (X)|, |S2 (X)|, |S3 (X)|} ≥ k and
P(E) ≥ 1 − 6k(0.86)2.5k .
For any fixed realization G ∼ H0 with set of clique vertices C = C1 ∪ C2 ∪ C3 (with corresponding supports Sk := Sk (X), k = 1, 2, 3), we generate a Gaussian random tensor Ỹ ∈ Rp×p×p
with independent entries such that
Ỹ (a, b, c) ∼ N (µ, 1)
where S1 = C1 , S2 = C2 −

if (a, b, c) ∈ S1 × S2 × S3 ;
p
2

Ỹ (a, b, c) ∼ N (0, 1)

and S3 = C3 − p. By Lemma 9, we have




2
dTV L Y (a, b, c) C , L Ỹ (a, b, c) C ≤ e(1−M )/2 ,
As a result, since M =

√

otherwise,

∀ a, b, c ∈ {1, 2, . . . , p}.

8 log N and p = N/3,



dTV L(Y), L(Ỹ) =EC dTV L(Y|C), L(Ỹ|C)
=EC

p
X

√





dTV L Y (a, b, c) C , L Ỹ (a, b, c) C



3 (1−M 2 )/2

≤p e

a,b,c=1

Now we show that L(Ỹ|E) is a mixture over



PX , X ∈ M0 (p, k, r, λ)

with λ =

Indeed, for any fixed C, let X̃ = E(Ỹ|C). Then,


X̃(a, b, c) = E Ỹ (a, b, c)|C = µ,

∀(a, b, c) ∈ S1 × S2 × S3 .

Recall that on E, min{|S1 |, |S2 |, |S3 |} ≥ k. Therefore, X̃ is of rank 1 and on E,




min σmin M1 (X̃) , σmin M2 (X̃) , σmin M3 (X̃)
24

≤

e
.
27N
3(1−τ )/4

p√
.
2 8 log 3p

≥µ
since µ =

1
2M .

p
bp3(1−τ )/4 c
|S1 ||S2 ||S3 | ≥ µk 3/2 ≥ µbp3(1−τ )/4 c = √
2 8 log 3p

The above fact indicates that X̃ ∈ M0 (p, k, r, λ). In other words, under H0 , for

any C conditioned on E, there exists X(C) ∈ M0 (p, k, r, λ) such that L(Ỹ|C) = PX(C) . Define

the probability distribution π0 = L X(C)|E supported on M0 (p, k, r, λ). Then L(Ỹ|E) = Pπ0
and



dTV L(Y), Pπ0 ≤dTV L(Y), L(Ỹ) + TV L(Ỹ), Pπ0
√
√
2.5k
e
e
c
≤
+ P(E ) ≤
+ 6k 0.86
.
27N
27N

C.9

Proof of Lemma 11

p
N/2 and N=3p. Let N1 := N/2 =

p
{1, 2, . . . , 3p
2 } and C1 ⊂ 1, 2, . . . , 2 , we have

Recall that κ ≤

3p
2 .

Since C is uniformly chosen from

Pκ/8 p/2 p 
p/2  p 

κ
κ
+
1
κ/8 7κ/8
s=0 s
P |C1 | ≤
≤
 κ−s ≤

N1
N1
8
8
κ
κ


κ + 1 κ (p)(p − 1) . . . (p − 7κ/8 + 1)(p/2)(p/2 − 1) . . . (p/2 − κ/8 + 1)
=
8
κ/8
(3p/2)(3p/2 − 1) . . . (3p/2 − κ + 1)
κ/8 κ + 1

κ+1
2 κ κ
≤
8e · 28 /38
≤
(8e)κ/8
=
(0.86)κ/8
8
3
8
8
 p 

n
k
where we used the fact p/2
s
κ−s increases for 0 ≤ s ≤ κ/8 and inequality k ≤ (ne/k) .
Therefore, with probability at least 1 −

κ+1
κ/8 ,
4 (0.86)

κ
7κ
≤ |C1 |, |C2 |, |C3 | ≤
.
8
8
Recall that κ = 20k and we conclude the proof.

25

