# Nonconvex Low-Rank Tensor Completion
from Noisy Data00footnotetext: Corresponding author: Yuxin Chen. A short
version of this work has appeared in NeurIPS 2019 [CLPC19].

Changxiao Cai Department of Electrical Engineering, Princeton University,
Princeton, NJ 08544, USA; Email<EMAIL_ADDRESS>Gen Li
Department of Electronic Engineering, Tsinghua University, Beijing, 10084,
China; Email<EMAIL_ADDRESS>H. Vincent Poor11footnotemark: 1
Yuxin Chen11footnotemark: 1

###### Abstract

We study a noisy tensor completion problem of broad practical interest,
namely, the reconstruction of a low-rank tensor from highly incomplete and
randomly corrupted observations of its entries. While a variety of prior work
has been dedicated to this problem, prior algorithms either are
computationally too expensive for large-scale applications, or come with sub-
optimal statistical guarantees. Focusing on “incoherent” and well-conditioned
tensors of a constant CP rank, we propose a two-stage nonconvex algorithm —
(vanilla) gradient descent following a rough initialization — that achieves
the best of both worlds. Specifically, the proposed nonconvex algorithm
faithfully completes the tensor and retrieves all individual tensor factors
within nearly linear time, while at the same time enjoying near-optimal
statistical guarantees (i.e. minimal sample complexity and optimal estimation
accuracy). The estimation errors are evenly spread out across all entries,
thus achieving optimal $\ell_{\infty}$ statistical accuracy. We have also
discussed how to extend our approach to accommodate asymmetric tensors. The
insight conveyed through our analysis of nonconvex optimization might have
implications for other tensor estimation problems.

Keywords: tensor completion, nonconvex optimization, gradient descent,
spectral methods, entrywise statistical guarantees, minimaxity

###### Contents

  1. 1 Introduction and motivation
    1. 1.1 Tensor completion from noisy entries
    2. 1.2 Computational and statistical challenges
  2. 2 Algorithm and main results
    1. 2.1 A two-stage nonconvex algorithm
    2. 2.2 Main results
    3. 2.3 Numerical experiments
    4. 2.4 Notation
  3. 3 Initialization
    1. 3.1 Step 1: subspace estimation via a spectral method
    2. 3.2 Step 2: retrieval of low-rank tensor factors from the subspace estimate
      1. 3.2.1 Procedure
      2. 3.2.2 Intuition
    3. 3.3 Other alternatives?
  4. 4 Related work
  5. 5 Analysis
    1. 5.1 Analysis for local convergence of GD
      1. 5.1.1 Preliminaries: gradient and Hessian calculation
      2. 5.1.2 Local strong convexity and smoothness
      3. 5.1.3 Leave-one-out gradient descent sequences
      4. 5.1.4 Key lemmas
    2. 5.2 Analysis for initialization: Part 1 (subspace estimation)
      1. 5.2.1 Key results
      2. 5.2.2 Leave-one-out sequences for subspace estimation
    3. 5.3 Analysis for initialization: Part 2 (retrieval of individual tensor factors)
      1. 5.3.1 Main results and leave-one-out sequences
      2. 5.3.2 Analysis
  6. 6 Discussion
  7. A Proofs for local convergence of GD
    1. A.1 Proof of Lemma 5.1
      1. A.1.1 Bounding $\alpha_{4}$
      2. A.1.2 Bounding $\alpha_{1}$
      3. A.1.3 Bounding $\alpha_{2}$
      4. A.1.4 Bounding $\alpha_{3}$
      5. A.1.5 Putting all this together
    2. A.2 Proof of Lemma 5.2
    3. A.3 Proof of Lemma 5.3
    4. A.4 Proof of Lemma 5.4
    5. A.5 Proof of Lemma 5.5
    6. A.6 Proof of Lemma 5.6
  8. B Proofs for retrieving tensor components
    1. B.1 Proof of Lemma 5.12
    2. B.2 Proof of Lemma 5.13
    3. B.3 Proof of Lemma B.1
      1. B.3.1 Controlling $\bm{F}^{\tau}$
      2. B.3.2 Controlling $\bm{C}^{\tau}$
    4. B.4 Proof of Lemma 5.14
      1. B.4.1 Controlling $\alpha_{3}$
      2. B.4.2 Controlling $\alpha_{2}$
      3. B.4.3 Controlling $\alpha_{1}$
      4. B.4.4 Combining $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$
    5. B.5 Proof of Lemma B.2
      1. B.5.1 Proximity of $\bm{M}^{\tau,\left(m\right)}$ and $\widehat{\bm{M}}^{\tau,\left(m\right)}$
      2. B.5.2 Proximity of $\bm{M}^{\tau}$ and $\widehat{\bm{M}}^{\tau,\left(m\right)}$
    6. B.6 Proof of Lemma 5.15
    7. B.7 Proof of Lemma B.3
    8. B.8 Proof of Lemma 5.16
      1. B.8.1 Controlling $\beta_{1}$
      2. B.8.2 Controlling $\beta_{2}$
      3. B.8.3 Controlling $\beta_{3}$
      4. B.8.4 Combining $\beta_{1}$, $\beta_{2}$ and $\beta_{3}$
    9. B.9 Proof of Lemma 5.17
    10. B.10 Proof of Corollary 5.11
  9. C Proof of Corollary 2.9
  10. D Auxiliary lemmas
    1. D.1 Statements of auxiliary lemmas
    2. D.2 Proof of Lemma D.1
    3. D.3 Proof of Lemma D.2 and Corollary D.3
      1. D.3.1 Proof of Lemma D.2
      2. D.3.2 Proof of Corollary D.3
    4. D.4 Proof of Lemma D.4
    5. D.5 Proof of Lemma D.5
    6. D.6 Proof of Lemma D.6
    7. D.7 Proof of Lemma D.7
    8. D.8 Proof of Lemma D.8
    9. D.9 Proof of Lemma D.9
    10. D.10 Proof of Lemma D.10
  11. E Extension to asymmetric tensors
    1. E.1 Problem settings
    2. E.2 Algorithms
    3. E.3 Numerical experiments
    4. E.4 Analysis ideas
    5. E.5 Proof of Lemma E.2

##  1 Introduction and motivation

###  1.1 Tensor completion from noisy entries

Estimation of low-complexity models from highly incomplete observations is a
fundamental task that spans a diverse array of science and engineering
applications. Arguably one of the most extensively studied problems of this
kind is matrix completion, where one wishes to recover a low-rank matrix given
only partial entries [DR16, CC18]. Moving beyond matrix-type data, a natural
higher-order generalization is low-rank tensor completion, which aims to
reconstruct a low-rank tensor when the vast majority of its entries are
unseen. There is certainly no shortage of applications that motivate the
investigation of tensor completion (e.g. personalized medicine [SN16, Paw19],
medical imaging [GRY11, SHKM14, CZA+17], seismic data analysis [KSS13,
EAHK13], multi-dimensional harmonic retrieval [CC14, YLW+17]). One concrete
example in operations research arises when learning the preference of
individual customers for a collection of products on the basis of historical
transactions [FL19, MP20]. Given the limited availability of transaction data
(e.g. each customer might only have purchased very few products before), it is
crucial to exploit multi-way customer-product interactions (e.g. users’
browsing and searching histories) in order to better predict the likelihood of
a customer purchasing a new product. Clearly, the presence of missing data and
the need of exploiting multi-way structure result in the task of tensor
completion. Additionally, tensor completion finds important applications in
visual data in-painting [LMWY13, LYX17], where one wishes to reconstruct video
data (or a sequence of images) from incomplete measurements. The video data
consist of at least two spatial variables and one temporal variable, whose
intrinsic connections are often modeled via certain low-complexity tensors.

For the sake of clarity, we phrase the problem formally before we proceed,
focusing on a simple model that already captures the intrinsic difficulty of
tensor completion in many aspects.111We focus on symmetric order-3 tensors
primarily for simplicity of presentation. Many of our findings naturally
extend to the more general case with asymmetric tensors of possibly higher
order. Detailed discussions are deferred to Appendix E due to the space
limits.  Imagine we are asked to estimate a symmetric order-three
tensor222Here, a tensor $\bm{T}\in\mathbb{R}^{d\times d\times d}$ is said to
be symmetric if $T_{j,k,l}=T_{k,j,l}=T_{k,l,j}=T_{l,k,j}=T_{j,l,k}=T_{l,j,k}$
for all $1\leq j,k,l\leq d$. $\bm{T}^{\star}\in\mathbb{R}^{d\times d\times d}$
from a small number of noisy entries

$\displaystyle
T_{j,k,l}=T_{j,k,l}^{\star}+E_{j,k,l},\qquad\forall(j,k,l)\in\Omega,$ (1)

where $T_{j,k,l}$ is the observed noisy entry at location $(j,k,l)$,
$E_{j,k,l}$ stands for the associated noise, and
$\Omega\subseteq\\{1,\cdots,d\\}^{3}$ is a symmetric index subset to sample
from. For notational simplicity, we set $\bm{T}=[T_{j,k,l}]_{1\leq j,k,l\leq
d}$ and $\bm{E}=[E_{j,k,l}]_{1\leq j,k,l\leq d}$, with $T_{j,k,l}=E_{j,k,l}=0$
for any $(j,k,l)\notin\Omega$. We adopt a random sampling model such that each
index $(j,k,l)$ ($j\leq k\leq l$) is included in $\Omega$ independently with
probability $p$. In addition, we know a priori that the unknown tensor
$\bm{T}^{\star}\in\mathbb{R}^{d\times d\times d}$ is a superposition of $r$
rank-one tensors (often termed canonical polyadic (CP) decomposition if $r$ is
minimal)

$\displaystyle\bm{T}^{\star}=\sum_{i=1}^{r}\bm{u}_{i}^{\star}\otimes\bm{u}_{i}^{\star}\otimes\bm{u}_{i}^{\star},\qquad\text{or
more concisely,}\qquad\bm{T}^{\star}=\sum_{i=1}^{r}\bm{u}_{i}^{\star\,\otimes
3},$ (2)

where each $\bm{u}_{i}^{\star}\in\mathbb{R}^{d}$ represents one of the $r$
low-rank tensor components / factors. Here and throughout, for any vectors
$\bm{a},\bm{b},\bm{c}\in\mathbb{R}^{d}$, the tensor
$\bm{a}\otimes\bm{b}\otimes\bm{c}$ is a $d\times d\times d$ array whose
$(j,k,l)$-th entry is given by $a_{j}b_{k}c_{l}$. The primary question is
this: can we hope to faithfully estimate $\bm{T}^{\star}$, as well as the
individual tensor factors $\\{\bm{u}_{i}^{\star}\\}_{1\leq i\leq r}$, from the
partially revealed entries (1), assuming that $r$ is reasonably small?

###  1.2 Computational and statistical challenges

Even though tensor completion conceptually resembles matrix completion in
various ways, it is considerably more challenging than the matrix counterpart.
This is perhaps not surprising, given that a plethora of natural tensor
problems (e.g. computing the spectral norm, finding the best low-rank
approximation) are all notoriously hard [HL13]. As a notable example, while
matrix completion is often efficiently solvable under nearly minimal sample
complexity [CR09, Gro11], all polynomial-time algorithms developed so far for
tensor completion — even in the noise-free case — require a sample size at
least exceeding the order of $rd^{3/2}$, which is substantially larger than
the degrees of freedom (i.e. $rd$) underlying the model (2). In fact, it is
widely conjectured that there exists a large computational barrier away from
the information-theoretic sampling limits [BM16].

With this fundamental gap in mind, the current paper focuses on the regime (in
terms of the sample size) that enables reliable tensor completion in
polynomial time. A variety of algorithms have been proposed that enjoy some
sort of theoretical guarantees in (at least part of) this regime, including
but not limited to spectral methods [MS18, CLC+20], sum-of-squares hierarchy
[BM16, PS17], nonconvex algorithms [JO14, XY17], and also convex relaxation
(based on proper unfolding) [GRY11, HMGW15, RPP13, GQ14]. While these are all
polynomial-time algorithms, most of the computational complexities supported
by prior theory remain prohibitively high when dealing with large-scale tensor
data — a point that we shall elaborate on later. The only exception is the
unfolding-based spectral method, which, however, fails to achieve exact
recovery as the noise vanishes. This leads to a critical question:

  * Q1: Is there any linear-time algorithm that is guaranteed to work for low-rank tensor completion?

Going beyond such computational concerns, one might naturally wonder whether
it is also possible for a fast algorithm to achieve a nearly un-improvable
statistical accuracy in the presence of noise. Towards this end, intriguing
stability guarantees have been established for sum-of-squares hierarchy in the
noisy settings [BM16], although this paradigm is computationally expensive for
large-scale data. The recent work [XYZ17] came up with a two-stage algorithm
(i.e. a spectral method followed by tensor power iterations) for noisy tensor
completion. Its estimation accuracy, however, falls short of achieving exact
recovery in the absence of noise. This gives rise to another question of
fundamental importance:

  * Q2: Can we achieve near-optimal statistical accuracy without compromising computational efficiency?

In this paper, we aim to address the above two questions by developing a
nonconvex algorithm that achieves optimal computational efficiency and
statistical accuracy all at once.

##  2 Algorithm and main results

###  2.1 A two-stage nonconvex algorithm

To address the above-mentioned challenges, a first impulse is to resort to the
following least squares problem:

$\displaystyle\underset{\bm{u}_{1,}\cdots,\bm{u}_{r}\in\mathbb{R}^{d}}{\text{minimize}}\quad\sum_{j,k,l\in\Omega}\Big{(}\left[\sum\nolimits_{i=1}^{r}\bm{u}_{i}^{\otimes
3}\right]_{j,k,l}-T_{j,k,l}\Big{)}^{2},$ (3)

or more concisely (up to proper re-scaling),

$\displaystyle\underset{\bm{U}\in\mathbb{R}^{d\times r}}{\text{minimize}}\quad
f(\bm{U}):=\frac{1}{6p}\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{i=1}^{r}\bm{u}_{i}^{\otimes
3}-\bm{T}\Big{)}\Big{\|}^{2}_{\mathrm{F}}$ (4)

if we take
$\bm{U}:=\left[\bm{u}_{1},\dots,\bm{u}_{r}\right]\in\mathbb{R}^{d\times r}$.
Here, we denote by $\mathcal{P}_{\Omega}(\bm{T})$ the orthogonal projection of
any tensor $\bm{T}$ onto the subspace of tensors which vanish outside of the
index set $\Omega$. This optimization problem, however, is highly nonconvex
(which involves minimizing a degree-6 polynomial), thus resulting in
computational intractability in general.

Fortunately, not all nonconvex problems are as daunting to solve as they may
seem. For example, recent years have seen a flurry of activity in low-rank
matrix factorization via nonconvex optimization, which provably achieves
optimal statistical accuracy and computational efficiency at once; see [CLC19]
for an overview of recent advances. Motivated by this strand of work, we
propose to solve (4) via a two-stage nonconvex paradigm, presented below in
reverse order. The whole procedure is summarized in Algorithms 1-3.

##### Gradient descent (GD).

Arguably one of the simplest optimization algorithms is gradient descent,
which adopts a gradient update rule

$\displaystyle\bm{U}^{t+1}=\bm{U}^{t}-\eta_{t}\nabla f(\bm{U}^{t}),\qquad
t=0,1,\cdots$ (5)

where $\eta_{t}$ is the learning rate or the stepsize, and
$\bm{U}^{t}\in\mathbb{R}^{d\times r}$ is the estimate in the $t$-th iteration.
The main computational burden in each iteration lies in gradient evaluation,
which, in this case, can be performed in time proportional to that taken to
read the data.

Despite the simplicity of this algorithm, two critical issues stand out and
might significantly affect its efficiency, which we shall bear in mind
throughout the algorithmic and theoretical development.

(i) Local stationary points and initialization. As is well known, GD is
guaranteed to find an approximate local stationary point, provided that the
learning rates do not exceed the inverse Lipschitz constant of the gradient
[Bub15]. There exist, however, local stationary points (e.g. saddle points or
spurious local minima) that might fall short of the desired statistical
properties. This requires us to properly avoid such undesired points, while
retaining computational efficiency. To address this issue, one strategy is to
first identify a rough initial guess within a local region surrounding the
global solution (which often helps rule out bad local minima), in order to
guarantee proper convergence of subsequent optimization procedures [LT17,
JO14]. As a side remark, while careful initialization might not be crucial for
several matrix recovery cases [CCFM19, GBW18, TV19], it does seem to be
critical in various tensor problems [RM14]. We shall elucidate this point in
Section 3.3.

(ii) Learning rates and regularization. Learning rates play a pivotal role in
determining the convergence properties of GD. The challenge, however, is that
the loss function (4) is overall not sufficiently smooth (i.e. its gradient
often has an exceedingly large Lipschitz constant), and hence generic
optimization theory recommends a pessimistically slow update rule (i.e. an
extremely small learning rate) so as to guard against over-shooting. This,
however, slows down the algorithm significantly, thus destroying the main
computational advantage of GD (i.e. low per-iteration cost). With this issue
in mind, prior literature suggests carefully designed regularization steps
(e.g. proper projection, regularized loss functions) in order to improve the
geometry of the optimization landscape [XY17]. In contrast, we argue that one
is allowed to take a constant learning rate — which is as aggressive as it can
possibly be — even without enforcing any regularization procedures.

Algorithm 1 Gradient descent for nonconvex tensor completion

1:Generate an initial estimate $\bm{U}^{0}\in\mathbb{R}^{d\times r}$ via
Algorithm 2.

2:for $t=0,1,\dots,t_{0}-1$ do

3: $\bm{U}^{t+1}=\bm{U}^{t}-\eta_{t}\nabla
f(\bm{U}^{t})=\bm{U}^{t}-\frac{\eta_{t}}{p}\mathcal{P}_{\Omega}\big{(}\sum\nolimits_{i=1}^{r}\big{(}\bm{u}_{i}^{t}\big{)}^{\otimes
3}-\bm{T}\big{)}\times_{1}^{\mathsf{seq}}\bm{U}^{t}\times_{2}^{\mathsf{seq}}\bm{U}^{t}$,
where $\times_{1}^{\mathsf{seq}}$ and $\times_{2}^{\mathsf{seq}}$ are defined
in Section 2.4.

##### Initialization.

Motivated by the above-mentioned issue (i), we develop a procedure that
guarantees a reasonable initial estimate. In a nutshell, the proposed
procedure consists of two steps:

  * (a)

Estimate the subspace spanned by the $r$ low-rank tensor factors
$\\{\bm{u}_{i}^{\star}\\}_{1\leq i\leq r}$ via a spectral method;

  * (b)

Disentangle individual low-rank tensor factors from this subspace estimate.

As we shall see momentarily, the total computational complexity of the
proposed initialization is $O(pd^{3})$ when $r=O(1)$, $\kappa=O(1)$ and $p\geq
1/d^{2}$ (where $\kappa$ is a sort of “condition number” defined later), which
is a linear-time algorithm. Note, however, that these two steps in the
initialization procedure are relatively more complicated to describe. To
improve the flow of the current paper, we postpone the details to Section 3.
The readers can catch a glimpse of these procedures in Algorithms 2-3.

Algorithm 2 Spectral initialization for nonconvex tensor completion

1:Let $\bm{U}\bm{\Lambda}\bm{U}^{\top}$ be the rank-$r$ eigen-decomposition of
$\displaystyle\bm{B}:=\mathcal{P}_{\mathsf{off}\text{-}\mathsf{diag}}(\bm{A}\bm{A}^{\top}),$
(6) where $\bm{A}=\mathsf{unfold}\big{(}p^{-1}\bm{T}\big{)}$ is the mode-1
matricization of $p^{-1}\bm{T}$, and
$\mathcal{P}_{\mathsf{off}\text{-}\mathsf{diag}}(\bm{Z})$ extracts out the
off-diagonal entries of $\bm{Z}$.

2:Output: an initial estimate $\bm{U}^{0}\in\mathbb{R}^{d\times r}$ on the
basis of $\bm{U}\in\mathbb{R}^{d\times r}$ using Algorithm 3.

Algorithm 3 Retrieval of low-rank tensor factors from a given subspace
estimate.

1:Input: number of restarts $L$, pruning threshold $\epsilon_{\mathsf{th}}$,
subspace estimate $\bm{U}\in\mathbb{R}^{d\times r}$ given by Algorithm 2.

2:for $\tau=1,\dots,L$ do

3: Generate an independent Gaussian vector
$\bm{g}^{\tau}\sim\mathcal{N}(0,\bm{I}_{d})$.

4:
$\big{(}\bm{\nu}^{\tau},\lambda_{\tau},\mathsf{spec}\text{-}\mathsf{gap}_{\tau}\big{)}\leftarrow\textsc{Retrieve-
one-tensor-factor}(\bm{T},p,\bm{U},\bm{g}^{\tau})$.

5:Generate tensor factor estimates
$\big{\\{}(\bm{w}^{1},\lambda_{1}),\dots,(\bm{w}^{r},\lambda_{r})\big{\\}}\leftarrow\textsc{Prune}(\big{\\{}\big{(}\bm{\nu}^{\tau},\lambda_{\tau},\mathsf{spec}\text{-}\mathsf{gap}_{\tau}\big{)}\big{\\}}_{\tau=1}^{L},\epsilon_{\mathsf{th}})$.

6:Output: initial estimate
$\bm{U}^{0}=\big{[}\lambda_{1}^{1/3}\bm{w}^{1},\dots,\lambda_{r}^{1/3}\bm{w}^{r}\big{]}.$

1:function Retrieve-one-tensor-factor($\bm{T},p,\bm{U},\bm{g}$)

2: Compute  $\displaystyle\bm{\theta}$
$\displaystyle=\bm{U}\bm{U}^{\top}\bm{g}=:\mathcal{P}_{\bm{U}}(\bm{g}),$ (7a)
$\displaystyle\bm{M}$ $\displaystyle=p^{-1}\bm{T}\times_{3}\bm{\theta},$ (7b)
where $\times_{3}$ is defined in Section 2.4.

3: Let $\bm{\nu}$ be the leading singular vector of $\bm{M}$ obeying
$\langle\bm{T},\bm{\nu}^{\otimes 3}\rangle\geq 0$, and set $\lambda=\langle
p^{-1}\bm{T},\bm{\nu}^{\otimes 3}\rangle$.

4: return
$\big{(}\bm{\nu},\lambda,\sigma_{1}(\bm{M})-\sigma_{2}(\bm{M})\big{)}$.

1:function
Prune($\big{\\{}\big{(}\bm{\nu}^{\tau},\lambda_{\tau},\mathsf{spec}\text{-}\mathsf{gap}_{\tau}\big{)}\big{\\}}_{\tau=1}^{L},\epsilon_{\mathsf{th}}$)

2: Set
$\Theta=\big{\\{}\big{(}\bm{\nu}^{\tau},\lambda_{\tau},\mathsf{spec}\text{-}\mathsf{gap}_{\tau}\big{)}\big{\\}}_{\tau=1}^{L}.$

3: for $i=1,\dots,r$ do

4: Choose
$(\bm{\nu}^{\tau},\lambda_{\tau},\mathsf{spec}\text{-}\mathsf{gap}_{\tau})$
from $\Theta$ with the largest $\mathsf{spec}\text{-}\mathsf{gap}_{\tau}$; set
$\bm{w}^{i}=\bm{\nu}^{\tau}$ and $\lambda_{i}=\lambda_{\tau}$.

5: Update
$\Theta\leftarrow\Theta\setminus\left\\{\big{(}\bm{\nu}^{\tau},\lambda_{\tau},\mathsf{spec}\text{-}\mathsf{gap}_{\tau}\big{)}\in\Theta:|\langle\bm{\nu}^{\tau},\bm{w}^{i}\rangle|>1-\epsilon_{\mathsf{th}}\right\\}$.

6: return
$\big{\\{}(\bm{w}^{1},\lambda_{1}),\dots,(\bm{w}^{r},\lambda_{r})\big{\\}}.$

###  2.2 Main results

Encouragingly, the proposed nonconvex algorithm provably achieves the best of
both worlds — in terms of statistical accuracy and computational efficiency —
for a class of low-rank, well-conditioned, and “incoherent” problem instances.
This subsection summarizes our main findings.

Before continuing, we note that one cannot hope to recover an arbitrary tensor
from highly sub-sampled and arbitrarily corrupted entries. In order to enable
provably valid recovery, the present paper focuses on a tractable model by
imposing the following assumptions.

######  Definition 2.1 (Incohrence and well-conditionedness).

Define the incoherence parameters and the condition number of $\bm{T}^{\star}$
as follows

$\displaystyle\mu_{0}$
$\displaystyle:=\frac{d^{3}\,\|\bm{T}^{\star}\|_{\infty}^{2}}{\|\bm{T}^{\star}\|_{\mathrm{F}}^{2}},$
(8a) $\displaystyle\mu_{1}$
$\displaystyle:=\frac{d\|\bm{u}_{i}^{\star}\|_{\infty}^{2}}{\,\|\bm{u}_{i}^{\star}\|_{2}^{2}},$
(8b) $\displaystyle\mu_{2}$
$\displaystyle:=\frac{d\langle\bm{u}_{i}^{\star},\bm{u}_{j}^{\star}\rangle^{2}}{\,\|\bm{u}_{i}^{\star}\|_{2}^{2}\|\bm{u}_{j}^{\star}\|_{2}^{2}},$
(8c) $\displaystyle\kappa$
$\displaystyle:=\frac{\max_{i}\|\bm{u}_{i}^{\star}\|_{2}}{\min_{i}\|\bm{u}_{i}^{\star}\|_{2}}.$
(8d)

######  Remark 2.2.

Here, $\mu_{0}$, $\mu_{1}$ and $\mu_{2}$ are termed the incoherence
parameters. Definitions (8a)-(8c) can be viewed as some sort of incoherence
conditions for the tensor. For instance, when $\mu_{0},\mu_{1}$ and $\mu_{2}$
are small, these conditions say that (1) the energy of tensor $\bm{T}^{\star}$
is (nearly) evenly spread across all entries; (2) each factor
$\bm{u}_{i}^{\star}$ is de-localized; (3) the factors
$\\{\bm{u}_{i}^{\star}\\}$ are nearly orthogonal to each other. Definition
(8d) is concerned with the “well-conditionedness” of the tensor, meaning that
each rank-1 component is of roughly the same size. In particular, we note that
an assumption on pairwise correlation (i.e. a constraint on $\mu_{2}$) is
often assumed in the literature of tensor decomposition / factorization (e.g.
[AGJ14, SLLC17, HZC20]).

For notational simplicity, we shall set

$\displaystyle\mu:=\max\left\\{\mu_{0},\,\mu_{1},\,\mu_{2}\right\\}.$ (9)

Note that our theory allows $\mu$ to grow with the problem dimension $d$ (in
fact, $\mu$ can be as large as $d/\mathrm{poly}\log(d)$).

######  Assumption 2.3 (Random noise).

Suppose that $\bm{E}$ is a symmetric random tensor, where
$\\{E_{j,k,l}\\}_{1\leq j\leq k\leq l\leq d}$ (cf. (1)) are independently
generated sub-Gaussian random variables with mean zero and variance
$\mathsf{Var}(E_{j,k,l})\leq\sigma^{2}$.

In addition, recognizing that there is a global permutational ambiguity issue
(namely, one cannot distinguish $\bm{u}_{1}^{\star},\cdots,\bm{u}_{r}^{\star}$
from an arbitrary permutation of them), we introduce the following loss
metrics to account for this ambiguity:

$\displaystyle\mathsf{dist}_{\mathrm{F}}(\bm{U},\bm{U}^{\star})$
$\displaystyle:=\min_{\bm{\Pi}\in\mathsf{perm}_{r}}\|\bm{U}\bm{\Pi}-\bm{U}^{\star}\|_{\mathrm{F}},$
(10a) $\displaystyle\mathsf{dist}_{\infty}(\bm{U},\bm{U}^{\star})$
$\displaystyle:=\min_{\bm{\Pi}\in\mathsf{perm}_{r}}\|\bm{U}\bm{\Pi}-\bm{U}^{\star}\|_{\infty},$
(10b) $\displaystyle\mathsf{dist}_{2,\infty}(\bm{U},\bm{U}^{\star})$
$\displaystyle:=\min_{\bm{\Pi}\in\mathsf{perm}_{r}}\|\bm{U}\bm{\Pi}-\bm{U}^{\star}\|_{2,\infty},$
(10c)

where $\mathsf{perm}_{r}$ stands for the set of $r\times r$ permutation
matrices. For notational simplicity, we also take

$\displaystyle\lambda_{\min}^{\star}:=\min_{1\leq i\leq
r}\left\|\bm{u}_{i}^{\star}\right\|_{2}^{3}\qquad\text{and}\qquad\lambda_{\max}^{\star}:=\max_{1\leq
i\leq r}\left\|\bm{u}_{i}^{\star}\right\|_{2}^{3}.$ (11)

With these notations in place, we are ready to present our main results. For
simplicity of presentation, we shall start with the setting where
$r,\mu,\kappa\asymp 1$.

######  Theorem 2.4.

Fix an arbitrary small constant $\delta>0$. Suppose that $r,\kappa,\mu=O(1)$,

$\displaystyle p\geq
c_{0}\frac{\log^{4}d}{d^{3/2}},\qquad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{d^{3/4}\log^{2}d},$ $\displaystyle
L=c_{2}\quad\text{and}\quad\epsilon_{\mathsf{th}}=c_{3}\bigg{(}\frac{\log
d}{d\sqrt{p}}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log^{2}d}{p}}+\sqrt{\frac{\log
d}{d}}\bigg{)}$

for some sufficiently large constants $c_{0},c_{2}>0$ and some sufficiently
small constants $c_{1},c_{3}>0$. The learning rate $\eta_{t}\equiv\eta$ is
taken to be a constant obeying $0<\eta\leq\lambda_{\min}^{\star
4/3}/\big{(}32\lambda_{\max}^{\star 8/3}\big{)}$. Then with probability at
least $1-\delta$,

$\displaystyle\mathsf{dist}_{\mathrm{F}}(\bm{U}^{t},\bm{U}^{\star})$
$\displaystyle\leq\Bigg{(}C_{1}\rho^{t}+C_{2}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\Bigg{)}\left\|\bm{U}^{\star}\right\|_{\mathrm{F}},$ (12a)
$\displaystyle\mathsf{dist}_{\infty}(\bm{U}^{t},\bm{U}^{\star})$
$\displaystyle\leq\mathsf{dist}_{2,\infty}(\bm{U}^{t},\bm{U}^{\star})\leq\Bigg{(}C_{3}\rho^{t}+C_{4}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\Bigg{)}\left\|\bm{U}^{\star}\right\|_{2,\infty}$ (12b)

hold simultaneously for all $0\leq t\leq t_{0}=d^{5}$. Here,
$0<C_{1},C_{3},\rho<1$ and $C_{2},C_{4}>0$ are some absolute constants.

######  Remark 2.5.

The theorem holds unchanged if $d^{5}$ is replaced by $d^{c}$ for an
arbitrarily large constant $c>0$.

######  Remark 2.6.

The upper bound $t_{0}$ on the iteration count arises from the leave-one-out
analysis when handling noisy observations. In short, the leave-one-out
argument can only provide high-probability bounds for each iteration, thus
requiring an upper bound on the iteration count if we desire a uniform bound
across iterations. Note that in the noiseless case, our results and analysis
hold for an arbitrarily large number of iterations.

As an immediate consequence of Theorem 2.4, we obtain appealing
$\ell_{\infty}$ statistical guarantees for estimating tensor entries, which
are previously rarely available (see Table 1). Specifically, let our tensor
estimate in the $t$-th iteration be

$\displaystyle\bm{T}^{t}:=\sum_{i=1}^{r}\bm{u}_{i}^{t}\otimes\bm{u}_{i}^{t}\otimes\bm{u}_{i}^{t},\qquad\text{where
}\bm{U}^{t}=[\bm{u}_{1}^{t},\cdots,\bm{u}_{r}^{t}]\in\mathbb{R}^{d\times r}.$
(13)

Then our result is this:

######  Corollary 2.7.

Fix an arbitrarily small constant $\delta>0$. Instate the assumptions of
Theorem 2.4. Then with probability at least $1-\delta$,

$\displaystyle\left\|\bm{T}^{t}-\bm{T}^{\star}\right\|_{\mathrm{F}}$
$\displaystyle\lesssim\left(C_{1}\rho^{t}+C_{2}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{T}^{\star}\right\|_{\mathrm{F}},$ (14a)
$\displaystyle\left\|\bm{T}^{t}-\bm{T}^{\star}\right\|_{\infty}$
$\displaystyle\lesssim\left(C_{3}\rho^{t}+C_{4}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{T}^{\star}\right\|_{\infty}$ (14b)

hold simultaneously for all $0\leq t\leq t_{0}=d^{5}$. Here,
$0<C_{1},C_{3},\rho<1$ and $C_{2},C_{4}>0$ are some absolute constants.

Several important implications are provided as follows. The discussion below
assumes $\lambda_{\max}^{\star}\asymp\lambda_{\min}^{\star}\asymp 1$ for
notational simplicity.

  1. 1.

Linear convergence. In the absence of noise, the proposed algorithm converges
linearly, namely, it provably attains $\varepsilon$ accuracy within
$O(\log(1/\varepsilon))$ iterations. Given the inexpensiveness of each
gradient iteration, this algorithm can be viewed as a linear-time algorithm,
which can almost be implemented as long as we can read the data. In the noisy
setting, the algorithm reaches an appealing statistical accuracy within a
logarithmic number of iterations.

  2. 2.

Near-optimal sample complexity. The fast convergence is guaranteed as soon as
the sample size exceeds the order of $d^{3/2}\mathrm{poly}\log d$. This
matches the minimal sample complexity — modulo some logarithmic factor — known
so far for any polynomial-time algorithm.

  3. 3.

Near-optimal statistical accuracy. The proposed algorithm converges
geometrically fast to a point with Euclidean error $O\big{(}\sigma\sqrt{(d\log
d)/p}\big{)}$. This matches the lower bound established in [XYZ17, Theorem 5]
up to some logarithmic factor, thus justifying the statistical optimality of
the proposed nonconvex algorithm.

  4. 4.

Entrywise estimation accuracy. In addition to the Euclidean statistical
guarantees, we have also established an entrywise error bound, which, to the
best of our knowledge, has not been established in any of the prior work. When
$t$ is sufficiently large, the iterates reach an entrywise error bound
$O\big{(}\sigma\sqrt{(\log d)/p}\big{)}$. This entrywise error bound is about
an order of $\sqrt{d}$ times smaller than the above $\ell_{2}$ error bound,
thereby implying that the estimation errors are evenly spread out across all
entries.

  5. 5.

Noise size. The above theory operates in the regime where
$\sigma\lesssim\sqrt{\frac{p}{d^{3/2}}}$ (modulo some log factor). Given that
we have $\|\bm{T}^{\star}\|_{\infty}\asymp d^{-3/2}$ in this case, our noise
size constraint can be equivalently written as (up to some log factor)

$\displaystyle\frac{\sigma}{\|\bm{T}^{\star}\|_{\infty}}\lesssim\sqrt{pd^{3/2}}.$
(15)

Since the sampling rate needs to satisfy $p\gg d^{-3/2}$, this condition
essentially allows the typical size of each noise component to be considerably
larger than the size of the corresponding entry of the truth, which covers a
broad range of practical scenarios.

  6. 6.

Implicit regularization. One appealing feature of our finding is the
simplicity of the iterative refinement stage of the algorithm. All of the
above statistical and computational benefits hold for vanilla gradient descent
(when properly initialized). This should be contrasted with prior work (e.g.
[XY17]) that relies on extra regularization terms to stabilize the
optimization landscape. In principle, vanilla gradient descent implicitly
constrains itself within a region of well-conditioned landscape, thus enabling
fast convergence without explicit regularization.

  7. 7.

No need of sample splitting. The theory developed herein does not require
fresh samples in each iteration. We note that sample splitting has been
frequently adopted in other context primarily to simplify mathematical
analysis. Nevertheless, it typically does not exploit the data in an efficient
manner (i.e. each data sample is used only once), thus resulting in the need
of a much larger sample size in practice.

| algorithm |  sample complexity |  computational complexity |  |  $\ell_{2}$ error
---
(noisy)
|  $\ell_{\infty}$ error
---
(noisy)
| recovery type
---
(noiseless)
| our
---
theory
spectral method \+ (vanilla) GD | $d^{1.5}$ | $pd^{3}$ | $\sigma\sqrt{\frac{d}{p}}$ | $\sigma\sqrt{\frac{1}{p}}$ | exact
[XYZ17] |  spectral initialization \+ tensor power method | $d^{1.5}$ | $pd^{3}$ | $(\|\bm{T}^{\star}\|_{\infty}+\sigma)\sqrt{\frac{d}{p}}$ | n/a | approximate
[XY17] |  spectral method \+ GD on manifold | $d^{1.5}$ | $\mathrm{poly}(d)$ | n/a | n/a | exact
[MS18] | spectral method | $d^{1.5}$ | $d^{3}$ | n/a | n/a | approximate
[BM16] | sum-of-squares | $d^{1.5}$ | $d^{15}$ | $\frac{\|\bm{T}^{\star}\|_{\mathrm{F}}}{\sqrt{pd^{1.5}}}+\sigma d^{1.5}$ | n/a | approximate
[PS17] | sum-of-squares | $d^{1.5}$ | $d^{10}$ | n/a | n/a | exact
| [YZ16]
---
[YZ17]
tensor nuclear norm minimization | $d$ | NP-hard | n/a | n/a | exact
Table 1: Comparison with prior theory for existing methods when
$r,\mu,\kappa\asymp 1$ (neglecting logarithmic factors).

We shall take a moment to discuss the merits of our approach in comparison to
prior work. One of the best-known polynomial-time algorithms is the degree-6
level of the sum-of-squares (SoS) hierarchy, which seems to match the
computationally feasible limit in terms of the sample complexity [BM16].
However, this approach has a well-documented limitation in that it involves
solving a semidefinite program of dimensions $d^{3}\times d^{3}$, which
requires enormous storage and computation power. The work [MS18] alleviates
this computational burden by resorting to a clever unfolding-based spectral
algorithm; it is a nearly linear-time procedure that enables near-minimal
sample complexity (among polynomial-time algorithms), although it does not
achieve exact recovery even in the absence of noise. The two-stage algorithm
developed by [XYZ17] — which is based on spectral initialization followed by
tensor power methods — shares similar advantages and drawbacks as [MS18].
Further, the recent work [XY17] proposes a polynomial-time nonconvex algorithm
based on gradient descent over Grassmann manifold (with a properly regularized
objective function), which is an extension of the nonconvex matrix completion
algorithm proposed by [KMO10a, KMO10b] to tensor data. The theory provided in
[XY17], however, does not provide explicit computational complexities. The
recent work [SY19] attempts tensor estimation via an interesting algorithm
adapted from collaborative filtering and investigates both $\ell_{2}$ and
$\ell_{\infty}$ estimation accuracy. This approach, however, does not
guarantee exact recovery in the absence of noise. We summarize and compare
several prior results in Table 1 (omitting logarithmic factors).

Thus far, we have concentrated on the low-rank, well-conditioned, and
incoherent case. Our main theory can be extended to cover a broader class of
scenarios, as stated below.

######  Theorem 2.8.

Fix an arbitrary small constant $\delta>0$. Suppose that $\kappa\asymp 1$,

$\displaystyle p\geq
c_{0}\frac{\mu^{4}r^{4}\log^{4}d}{d^{3/2}},\qquad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{\mu r^{3/2}d^{3/4}\log^{2}d},\quad r\leq
c_{2}\left(\frac{d}{\mu^{6}\log^{6}d}\right)^{1/6},$ $\displaystyle
L=c_{3}r^{2\kappa^{2}}\log^{3/2}r\qquad\text{and}\qquad\epsilon_{\mathsf{th}}=c_{4}\bigg{(}\frac{\mu
r\log
d}{d\sqrt{p}}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{rd\log^{2}d}{p}}+\sqrt{\frac{\mu
r\log d}{d}}\bigg{)}$

for some sufficiently large constants $c_{0},c_{3}>0$ and some sufficiently
small constants $c_{1},c_{2},c_{4}>0$. The learning rate $\eta_{t}\equiv\eta$
is taken to be a constant obeying $0<\eta\leq\lambda_{\min}^{\star
4/3}/\big{(}32\lambda_{\max}^{\star 8/3}\big{)}$. Then with probability at
least $1-\delta$,

$\displaystyle\mathsf{dist}_{\mathrm{F}}(\bm{U}^{t},\bm{U}^{\star})$
$\displaystyle\leq\left(C_{1}\rho^{t}+C_{2}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{\mathrm{F}}$ (16a)
$\displaystyle\mathsf{dist}_{\infty}(\bm{U}^{t},\bm{U}^{\star})$
$\displaystyle\leq\mathsf{dist}_{2,\infty}(\bm{U}^{t},\bm{U}^{\star})\leq\left(C_{3}\rho^{t}+C_{4}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{2,\infty}$ (16b)

hold simultaneously for all $0\leq t\leq t_{0}=d^{5}$. Here,
$0<C_{1},C_{3},\rho<1$ and $C_{2},C_{4}>0$ are some absolute constants.

######  Corollary 2.9.

Fix an arbitrarily small constant $\delta>0$. Instate the assumptions of
Theorem 2.8. Then with probability at least $1-\delta$,

$\displaystyle\left\|\bm{T}^{t}-\bm{T}^{\star}\right\|_{\mathrm{F}}$
$\displaystyle\lesssim\left(C_{1}\rho^{t}+C_{2}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{T}^{\star}\right\|_{\mathrm{F}},$ (17a)
$\displaystyle\left\|\bm{T}^{t}-\bm{T}^{\star}\right\|_{\infty}$
$\displaystyle\lesssim\left(C_{3}\rho^{t}+C_{4}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{\mu^{3}rd\log
d}{p}}\right)\left\|\bm{T}^{\star}\right\|_{\infty}$ (17b)

hold simultaneously for all $0\leq t\leq t_{0}=d^{5}$. Here,
$0<C_{1},C_{3},\rho<1$ and $C_{2},C_{4}>0$ are some absolute constants.

######  Remark 2.10.

Clearly, Theorem 2.8 and Corollary 2.9 subsume Theorem 2.4 and Corollary 2.7
as a special case respectively.

######  Remark 2.11.

Our theorems require the rank $r$ to not exceed $o(d^{1/6})$, which, we
believe, is an artifact of the current nonconvex analysis (particularly for
the initialization stage). For instance, our local convergence analysis is
built upon strong convexity and smoothness, which holds only within a
sufficiently small neighborhood surrounding the truth; given that the diameter
of this neighborhood is no more than $o(1/r)$, our analysis requires an
initial guess with higher accuracy than expected, thus leading to our rank
constraint. It might be possible to improve the rank dependency via more
refined analysis, and we leave it to future investigation.

In a nutshell, this theorem reveals intriguing theoretical support (including
both $\ell_{\mathrm{F}}$ and $\ell_{2,\infty}$ bounds) for more general
settings. Assuming that the condition number $\kappa\asymp 1$, the nonconvex
algorithm we propose is guaranteed to succeed in polynomial time. Note,
however, that our theoretical dependency (including both sample and
computational complexities) on the rank $r$ and the incoherence parameter
$\mu$ are likely loose and sub-optimal. In addition, if $\kappa$ is allowed to
grow with $d$, then the current theory requires a large number of restart
attempts during the initialization stage, resulting in a very high
computational burden. Improving these aspects, however, calls for a much more
refined analysis framework, which we leave for future investigation.

###  2.3 Numerical experiments

We carry out a series of numerical experiments to corroborate our theoretical
findings. Before proceeding, recall that Theorem 2.8 only guarantees
successful recovery with probability $1-\delta$ for some small constant
$\delta$; this means that we shall not anticipate a very high success rate
(e.g. $1-O(d^{-5})$) as in the matrix recovery case. As we shall make clear
shortly, this happens mainly because the initialization stage works only with
probability $1-\delta$, where the uncertainty largely depends on the random
vectors $\\{\bm{g}^{\tau}\\}_{1\leq\tau\leq L}$. With this observation in
mind, we recommend the following modification to improve the empirical success
rate:

  * •

Run Algorithm 2 independently for $t_{\mathsf{init}}=5$ times to obtain
multiple initial estimates (denoted by
$\bm{U}^{0}_{[1]},\cdots,\bm{U}^{0}_{[t_{\mathsf{init}}]}$); select the one
achieving the smallest empirical loss, namely

$\displaystyle\bm{U}^{0}_{\mathsf{best}}=\underset{{\bm{U}\in\\{\bm{U}^{0}_{[i]}\\}}_{1\leq
i\leq t_{\mathsf{init}}}}{\arg\min}f(\bm{U}).$ (18)

  * •

Run Algorithm 1 with the initial point $\bm{U}^{0}$ set to be
$\bm{U}^{0}_{\mathsf{best}}$.

The final estimates for the low-rank factor and the whole tensor are denoted
respectively by

$\displaystyle\widehat{\bm{U}}=\bm{U}^{t_{0}}\qquad\text{and}\qquad\widehat{\bm{T}}=\sum\nolimits_{i=1}^{r}\bm{u}^{t_{0}}_{i}\otimes\bm{u}^{t_{0}}_{i}\otimes\bm{u}^{t_{0}}_{i},$
(19)

where
$\bm{U}^{t_{0}}=[\bm{u}^{t_{0}}_{1},\cdots,\bm{u}^{t_{0}}_{r}]\in\mathbb{R}^{d\times
r}$ is the iterate returned by Algorithm 1, with $t_{0}$ the total number of
gradient iterations. In the sequel, we generate the true tensor
$\bm{T}^{\star}=\sum_{1\leq i\leq r}\bm{u}_{i}^{\star\,\otimes 3}$ randomly in
such a way that
$\bm{u}_{i}^{\star}\overset{\mathrm{i.i.d.}}{\sim}\mathcal{N}(\bm{0},\bm{I}_{d})$.
The learning rates are taken to be $\eta_{t}\equiv 0.2$ unless otherwise
noted.

We start with numerical convergence rates of our algorithm in the absence of
noise. Set $d=100$, $r=4$, $p=0.1$, $L=16$ and $\epsilon_{\mathsf{th}}=0.4$.
Figure 1(a) the numerical estimation errors vs. iteration count $t$ in a
typical Monte Carlo trial. Here, four kinds of estimation errors are reported:
(1) the relative Frobenius norm error
$\frac{\mathsf{dist}_{\mathrm{F}}(\bm{U}^{t},\bm{U}^{\star})}{\|\bm{U}^{\star}\|_{\mathrm{F}}}$;
(2) the relative $\|\cdot\|_{2,\infty}$ error
$\frac{\mathsf{dist}_{2,\infty}(\bm{U}^{t},\bm{U}^{\star})}{\|\bm{U}^{\star}\|_{2,\infty}}$;
(3) the relative Frobenius norm error
$\frac{\|\bm{T}^{t}-\bm{T}^{\star}\|_{\mathrm{F}}}{\|\bm{T}^{\star}\|_{\mathrm{F}}}$;
(4) the relative $\ell_{\infty}$ error
$\frac{\|\bm{T}^{t}-\bm{T}^{\star}\|_{\infty}}{\|\bm{T}^{\star}\|_{\infty}}$.
Here,
$\bm{T}^{t}=\sum_{i=1}^{r}\bm{u}_{i}^{t}\otimes\bm{u}_{i}^{t}\otimes\bm{u}_{i}^{t}$
with $\bm{U}^{t}=[\bm{u}^{t}_{1},\cdots,\bm{u}^{t}_{r}]$. For all these
metrics, the numerical estimation errors decay geometrically fast.

Next, we study the phase transition (in terms of the success rates for exact
recovery) in the noise-free settings. Set $d=100$, $r=4$, $L=16$ and
$\epsilon_{\mathsf{th}}=0.4$. For the sake of comparisons, we also report the
numerical performance of the tensor power method (TPM) followed by gradient
descent. When running the tensor power method, we set both the number of
iterations and the restart number to be $16$. Each trial is claimed to succeed
if the relative $\ell_{2}$ error obeys
$\frac{\mathsf{dist}_{\mathrm{F}}(\widehat{\bm{U}},\bm{U}^{\star})}{\|\bm{U}^{\star}\|_{\mathrm{F}}}\leq
0.01$. Figure 1(b) plots the empirical success rates over 100 independent
Monte Carlo trials. As can be seen, our initialization algorithm outperforms
the tensor power method.

The third series of experiments is concerned with the dependence of the
success rate on the rank $r$. Let us set $p=rd^{-3/2}\log^{2}d$, $L=r^{2}$ and
$\epsilon_{\mathsf{th}}=0.4$, and the success recovery criterion is the same
as above. Figure 1(c) depicts the empirical success rates (over $100$
independent Monte Carlo trials) as the rank $r$ varies. As can be seen from
the plots, the proposed algorithm is able to achieve exact reconstruction as
long as the rank $r$ is sufficiently small compared to $d$. The plausible
range of $r$, however, seems and seems to be larger than our theoretic
requirement $r=o(d^{1/6})$. This, once again, suggests the need of future
investigation to pin down the best possible dependency on $r$.

Finally, we consider the numerical estimation accuracy of our algorithm. Take
$t_{0}=100$, $d=100$, $r=4$, $p=0.1$, $L=16$ and $\epsilon_{\mathsf{th}}=0.4$.
Define the signal-to-noise ratio (SNR) to be
$\mathrm{SNR}=\frac{\left\|\bm{T}^{\star}\right\|_{\mathrm{F}}^{2}/d^{3}}{\sigma^{2}}$.
We report in Figure 1(d) three types of squared relative errors (namely,
$\frac{\mathsf{dist}^{2}_{\mathrm{F}}(\widehat{\bm{U}},\bm{U}^{\star})}{\|\bm{U}^{\star}\|^{2}_{\mathrm{F}}}$,
$\frac{\mathsf{dist}^{2}_{2,\infty}(\widehat{\bm{U}},\bm{U}^{\star})}{\|\bm{U}^{\star}\|^{2}_{2,\infty}}$
and
$\frac{\|\widehat{\bm{T}}-\bm{T}^{\star}\|^{2}_{\infty}}{\|\bm{T}^{\star}\|^{2}_{\infty}}$)
vs. SNR. Figure 1(d) illustrates that all three types of relative squared
errors scale inversely proportional to the SNR (since the slope in the figure
is roughly $-1$), which is consistent with our statistical guarantees.

|  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 1: (a) relative errors of the estimates $\bm{U}^{t}$ and $\bm{T}^{t}$
vs. iteration count $t$ for noiseless tensor completion, where $d=100$, $r=4$,
$p=0.1$; (b) empirical success rate vs. sampling rate, where $d=100$, $r=4$;
(c) empirical success rate vs. rank, where $p=rd^{-3/2}\log^{2}d$; (d) squared
relative errors vs. SNR for noisy settings, where $d=100$, $r=4$ and $p=0.1$.
Each point in (b), (c) and (d) is averaged over $100$ independent Monte Carlo
trials.

###  2.4 Notation

Before proceeding, we gather a few notations that will be used throughout this
paper. First of all, for any matrix $\bm{M}\in\mathbb{R}^{d\times d}$, we let
$\left\|\bm{M}\right\|$ and $\left\|\bm{M}\right\|_{\mathrm{F}}$ denote the
operator norm (or the spectral norm) and the Frobenius norm of $\bm{M}$,
respectively, and let $\bm{M}_{i,:}$ and $\bm{M}_{:,i}$ denote the $i$-th row
and $i$-th column, respectively. In addition,
$\lambda_{1}(\bm{M})\geq\lambda_{2}(\bm{M})\geq\dots\geq\lambda_{d}(\bm{M})$
denote the eigenvalues of $\bm{M}$ and
$\sigma_{1}(\bm{M})\geq\sigma_{2}(\bm{M})\geq\dots\geq\sigma_{d}(\bm{M})$
denote the singular values of $\bm{M}$.

For any tensor $\bm{T}\in\mathbb{R}^{d\times d\times d}$, let
$\bm{T}_{i,:,:}\in\mathbb{R}^{d\times d}$ denote the mode-$1$ $i$-slice with
entries $\left(\bm{T}_{i,:,:}\right)_{j,k}=T_{i,j,k}$, and $\bm{T}_{:,i,:}$
and $\bm{T}_{:,:,i}$ are defined in a similar way. For any tensors
$\bm{T},\bm{R}\in\mathbb{R}^{d\times d\times d}$, the inner product is defined
as $\langle\bm{T},\bm{R}\rangle:=\sum_{j,k,l}T_{j,k,l}R_{j,k,l}$. The
Frobenius norm of $\bm{T}$ is defined as
$\left\|\bm{T}\right\|_{\mathrm{F}}:=\sqrt{\langle\bm{T},\bm{T}\rangle}$. For
any vectors $\bm{u},\bm{v}\in\mathbb{R}^{d}$, we define the vector products of
a tensor $\bm{T}\in\mathbb{R}^{d\times d\times d}$ — denoted by
$\bm{T}\times_{3}\bm{u}\in\mathbb{R}^{d\times d}$ and
$\bm{T}\times_{1}\bm{u}\times_{2}\bm{v}\in\mathbb{R}^{d}$ — such that

$\displaystyle\big{[}\bm{T}\times_{3}\bm{u}\big{]}_{ij}\leavevmode\nobreak\ $
$\displaystyle:=\leavevmode\nobreak\ \sum\nolimits_{1\leq k\leq
d}T_{i,j,k}u_{k},\qquad$ $\displaystyle 1\leq i,j\leq d;$ (20a)
$\displaystyle\big{[}\bm{T}\times_{1}\bm{u}\times_{2}\bm{v}\big{]}_{k}\leavevmode\nobreak\
$ $\displaystyle:=\leavevmode\nobreak\ \sum\nolimits_{1\leq i,j\leq
d}T_{i,j,k}u_{i}v_{j},\qquad$ $\displaystyle 1\leq k\leq d.$ (20b)

The products $\bm{T}\times_{2}\bm{u}\in\mathbb{R}^{d\times d}$,
$\bm{T}\times_{3}\bm{u}\in\mathbb{R}^{d\times d}$,
$\bm{T}\times_{1}\bm{u}\times_{3}\bm{v}\in\mathbb{R}^{d}$ and
$\bm{T}\times_{2}\bm{u}\times_{3}\bm{v}\in\mathbb{R}^{d}$ are defined in a
similar manner. For any
$\bm{U}=[\bm{u}_{1},\cdots,\bm{u}_{r}]\in\mathbb{R}^{d\times r}$ and
$\bm{V}=[\bm{v}_{1},\cdots,\bm{v}_{r}]\in\mathbb{R}^{d\times r}$, we further
define

$\displaystyle\bm{T}\times_{1}^{\mathsf{seq}}\bm{U}\times_{2}^{\mathsf{seq}}\bm{V}:=[\bm{T}\times_{1}\bm{u}_{i}\times_{2}\bm{v}_{i}]_{1\leq
i\leq r}\in\mathbb{R}^{d\times r}.$ (21)

In addition, the operator norm of $\bm{T}$ is defined as

$\displaystyle\left\|\bm{T}\right\|:=\sup_{\bm{u},\bm{v},\bm{w}\,\in\,\mathbb{S}^{d-1}}\left\langle\bm{T},\bm{u}\otimes\bm{v}\otimes\bm{w}\right\rangle,$
(22)

where $\mathbb{S}^{d-1}:=\\{\bm{u}\in\mathbb{R}^{d}\mid\|\bm{u}\|_{2}=1\\}$
indicates the unit sphere in $\mathbb{R}^{d}$.

Further, $f(n)\lesssim g(n)$ or $f(n)=O(g(n))$ means that $|f(n)/g(n)|\leq
C_{1}$ for some constant $C_{1}>0$; $f(n)\gtrsim g(n)$ means that
$|f(n)/g(n)|\geq C_{2}$ for some constant $C_{2}>0$; $f(n)\asymp g(n)$ means
that $C_{1}\leq|f(n)/g(n)|\leq C_{2}$ for some constants $C_{1},C_{2}>0$;
$f(n)=o(g(n))$ means that $\lim_{n\rightarrow\infty}f(n)/g(n)=0$. In addition,
$f(n)\ll g(n)$ means that $f(n)\leq c_{1}g(n)$ for some sufficiently small
constant $c_{1}>0$, and $f(n)\gg g(n)$ means that $f(n)\geq c_{2}g(n)$ for
some sufficiently large constant $c_{2}>0$.

##  3 Initialization

This section presents formal details of the proposed two-step initialization,
accompanied by some intuition. Recall that the proposed initialization
procedure consists of two steps, which we discuss separately.

###  3.1 Step 1: subspace estimation via a spectral method

The spectral algorithm is often applied in conjunction with simple “unfolding”
(or “matricization”) to estimate the subspace spanned by the $r$ factors
$\\{\bm{u}_{i}^{\star}\\}_{1\leq i\leq r}$. This strategy is partly motivated
by prior approaches developed for covariance estimation with missing data
[Lou14, MS18, CLC+20]. We provide a brief introduction below.

Let

$\displaystyle\bm{A}=\mathsf{unfold}^{1\times
2}\big{(}\tfrac{1}{p}\bm{T}\big{)}\in\mathbb{R}^{d\times d^{2}},\quad\text{or
more concisely}\leavevmode\nobreak\ \leavevmode\nobreak\
\bm{A}=\mathsf{unfold}\big{(}\tfrac{1}{p}\bm{T}\big{)}\in\mathbb{R}^{d\times
d^{2}}$ (23)

be the mode-1 matricization of $p^{-1}\bm{T}$ (namely,
$\frac{1}{p}T_{i,j,k}=A_{i,(j-1)d+k}$ for any $1\leq i,j,k\leq d$) [KB09]. The
rationale of this step is that: under our model, the unfolded matrix $\bm{A}$
obeys

$\displaystyle\mathbb{E}[\bm{A}]=\mathsf{unfold}\big{(}\bm{T}^{\star}\big{)}=\sum_{i=1}^{r}\bm{u}_{i}^{\star}\left(\bm{u}_{i}^{\star}\otimes\bm{u}_{i}^{\star}\right)^{\top}=:\bm{A}^{\star},$
(24)

whose column space is precisely the span of $\\{\bm{u}^{\star}\\}_{1\leq i\leq
r}$. This motivates one to estimate the $r$-dimensional column space of
$\mathbb{E}[\bm{A}]$ from $\bm{A}$. Towards this, a natural strategy is to
look at the principal subspace of $\bm{A}\bm{A}^{\top}$. However, the diagonal
entries of $\bm{A}\bm{A}^{\top}$ bear too much influence on the principal
directions and need to be properly down-weighed. The current paper chooses to
work with the principal subspace of the following matrix that zeros out all
diagonal components:

$\displaystyle\bm{B}:=\mathcal{P}_{\mathsf{off}\text{-}\mathsf{diag}}(\bm{A}\bm{A}^{\top}),$
(25)

where $\mathcal{P}_{\mathsf{off}\text{-}\mathsf{diag}}(\bm{Z})$ extracts out
the off-diagonal entries of a squared matrix $\bm{Z}$. If we let
$\bm{U}\in\mathbb{R}^{d\times r}$ be an orthonormal matrix whose columns are
the top-$r$ eigenvectors of $\bm{B}$, then $\bm{U}$ serves as our subspace
estimate. See Algorithm 2 for a summary of the procedure.

###  3.2 Step 2: retrieval of low-rank tensor factors from the subspace
estimate

####  3.2.1 Procedure

As it turns out, it is possible to obtain rough (but reasonable) estimates of
all individual low-rank tensor factors $\\{\bm{u}_{i}^{\star}\\}_{1\leq i\leq
r}$ — up to global permutation — given a reliable subspace estimate $\bm{U}$.
This is in stark contrast to the low-rank matrix recovery case, where there
exists some global rotational ambiguity that prevents us from disentangling
the $r$ factors of interest.

We begin by describing how to retrieve one tensor factor from the subspace
estimate — a procedure summarized in $\textsc{Retrieve-one-tensor-factor}($).
Let us generate a random vector from the provided subspace $\bm{U}$ (which has
orthonormal columns), that is,

$\displaystyle\bm{\theta}=\underset{\text{projection of }\bm{g}\text{ onto
}\bm{U}}{\underbrace{\bm{U}\bm{U}^{\top}\bm{g}}},\qquad\bm{g}\sim\mathcal{N}(\bm{0},\bm{I}_{d}).$
(26)

The rescaled tensor data $p^{-1}\bm{T}$ is then transformed into a matrix via
proper “projection” along this random direction $\bm{\theta}$, namely,

$\displaystyle\bm{M}=\tfrac{1}{p}\bm{T}\times_{3}\bm{\theta}\in\mathbb{R}^{d\times
d}.$ (27)

Our estimate for a tensor factor is then given by $\lambda^{1/3}\bm{\nu}$,
where $\bm{\nu}$ is the leading singular vector of $\bm{M}$ obeying
$\langle\bm{T},\bm{\nu}^{\otimes 3}\rangle\geq 0$, and $\lambda$ is taken as
$\lambda=\big{\langle}p^{-1}\bm{T},\bm{\nu}^{\otimes 3}\big{\rangle}$.
Informally, $\bm{\nu}$ reflects the direction of the component
$\bm{u}_{i}^{\star}$ that exhibits the largest correlation with the random
direction $\bm{\theta}$, and $\lambda$ forms an estimate of the corresponding
size $\|\bm{u}^{\star}_{i}\|_{2}$.

A challenge remains, however, as there are oftentimes more than one tensor
factors to estimate. To address this issue, we propose to re-run the
aforementioned procedure multiple times, so as to ensure that we get to
retrieve each tensor factor of interest at least once. We will then apply a
careful pruning procedure (i.e. $\textsc{Prune}($)) to remove redundancy.

####  3.2.2 Intuition

To develop some intuition about the above procedure, consider the “heuristic”
case where
$\bm{\theta}=\bm{U}^{\star}(\bm{U}^{\star\top}\bm{U}^{\star})^{-1}\bm{U}^{\star\top}\bm{g}$,
namely, the idealistic scenario where the subspace estimate $\bm{U}$ is
accurate. Averaging out the randomness in the sampling pattern and the noise,
we see that the expected projected matrix (27) takes the following form:

$\displaystyle\mathbb{E}\big{[}\bm{M}\mid\bm{\theta}\big{]}=\bm{T}^{\star}\times_{3}\bm{\theta}=\sum_{i=1}^{r}\langle\bm{\theta},\bm{u}_{i}^{\star}\rangle\bm{u}_{i}^{\star}\bm{u}_{i}^{\star\top}.$

As a result, in the incoherent case where $\\{\bm{u}_{j}^{\star}\\}$ are
nearly orthogonal to each other, the leading singular vector of
$\mathbb{E}\big{[}\bm{M}\mid\bm{\theta}\big{]}$ — and hence that of $\bm{M}$
(i.e. $\bm{w}$) — is expected to be reasonably close to the factor
$\bm{u}_{i}^{\star}$ that enjoys the largest projected coefficient. In other
words, we expect

$\bm{\nu}\approx\frac{1}{\left\|\bm{u}_{i}^{\star}\right\|_{2}}\bm{u}_{i}^{\star},\qquad\text{where
}i=\arg\max_{1\leq j\leq
r}\big{|}\langle\bm{\theta},\bm{u}_{j}^{\star}\rangle\big{|}.$ (28)

In the mean time, armed with (28) and the incoherence assumption (such that
$\bm{u}_{i}^{\star}$ and $\bm{u}_{j}^{\star}$ are nearly orthogonal for $i\neq
j$), one might have

$\displaystyle\lambda=\big{\langle}\bm{T}^{\star},\bm{\nu}^{\otimes
3}\big{\rangle}\approx\frac{1}{\left\|\bm{u}_{i}^{\star}\right\|_{2}^{3}}\left\langle\bm{T}^{\star},\bm{u}_{i}^{\star\otimes
3}\right\rangle\approx\frac{1}{\left\|\bm{u}_{i}^{\star}\right\|_{2}^{3}}\left\langle\bm{u}_{i}^{\star\otimes
3},\bm{u}_{i}^{\star\otimes
3}\right\rangle=\left\|\bm{u}_{i}^{\star}\right\|_{2}^{3},$ (29)

thus explaining our choice of $\lambda$ in the proposed procedure. These
arguments hint at the ability of our procedure in retrieving one tensor factor
in each round.

The above intuitive argument, however, does not explain why we need to first
project a random vector $\bm{g}$ onto the (approximate) column space of
$\bm{U}^{\star}$. While we won’t go into detailed calculations here, we remark
in passing a crucial high variability issue: without proper projection, the
perturbation incurred by both the missing data and the noise might far exceed
the strength of the true signal. As a result, it is advised to first project
the data onto the desired subspace, in the hope of amplifying the signal-to-
noise ratio.

###  3.3 Other alternatives?

The careful reader may naturally wonder whether a careful initialization is
pivotal in achieving fast convergence. While a thorough answer to this has yet
to be developed, we shall point out some alternatives that seem sub-optimal in
both theory and practice. To simplify the presentation, the current subsection
focuses on the rank-1 noiseless case, where

$\displaystyle\bm{T}^{\star}=\bm{u}^{\star\,\otimes
3},\qquad\bm{T}=\tfrac{1}{p}\mathcal{P}_{\Omega}(\bm{T}^{\star}),\qquad\|\bm{u}^{\star}\|_{2}=1.$
(30)

Since the decision variable is now a $d$-dimensional vector, we shall employ
the conventional notation $\bm{u}^{t}$ to represent $\bm{U}^{t}$.

##### Random initialization.

We find it instrumental to begin with the population-level analysis, which
corresponds to the scenario with no missing data and noise ($p=1$ and
$\sigma=0$). A little calculation gives

$\displaystyle\mathbb{E}\left[\bm{u}^{1}\mid\bm{u}^{0}\right]$
$\displaystyle=\mathbb{E}\left[\bm{u}^{0}-\eta\nabla
f(\bm{u}^{0})\mid\bm{u}^{0}\right]=\left(1-\eta\|\bm{u}^{0}\|_{2}^{4}\right)\bm{u}^{0}+\eta\langle\bm{u}^{0},\bm{u}^{\star}\rangle^{2}\bm{u}^{\star}.$
(31)

As an immediate consequence, the expected correlation between the next iterate
and the truth obeys

$\displaystyle\mathbb{E}\left[\langle\bm{u}^{1},\bm{u}^{\star}\rangle\mid\bm{u}^{0}\right]$
$\displaystyle=\left\\{1-\eta\|\bm{u}^{0}\|_{2}^{4}+\eta\langle\bm{u}^{0},\bm{u}^{\star}\rangle\|\bm{u}^{\star}\|_{2}^{2}\right\\}\langle\bm{u}^{0},\bm{u}^{\star}\rangle.$

This means that if $\bm{u}^{0}$ and $\bm{u}^{\star}$ are positively correlated
and if the initial guess $\bm{u}^{0}$ is sufficiently small,333In fact, if a
random initialization $\bm{u}^{0}$ is not small, then one can easily show
that, with high probability, the $\ell_{2}$ norm of $\bm{u}^{t}$ is going to
drop geometrically fast at the beginning.  then one has

$\displaystyle\mathbb{E}\left[\langle\bm{u}^{1},\bm{u}^{\star}\rangle\mid\bm{u}^{0}\right]$
$\displaystyle\approx\left(1+\eta\langle\bm{u}^{0},\bm{u}^{\star}\rangle\|\bm{u}^{\star}\|_{2}^{2}\right)\langle\bm{u}^{0},\bm{u}^{\star}\rangle;$
(32)

a similar recursion holds for $\bm{u}^{t}$. As a result, the GD iterates are
expected to get increasingly more aligned with the truth, at least at the
population level. Caution needs to be exercised, however, that this
population-level analysis alone fails to capture what is happening in the
finite-sample case. In what follows, we point out potential issues with random
initialization.

Consider the case where $\bm{u}^{0}$ is generated as a vector of i.i.d.
Gaussian random variables. Suppose that $\bm{u}^{0}$ and $\bm{u}^{\star}$ are
positively correlated and that $\|\bm{u}^{0}\|_{2}$ is sufficiently small. It
is easily seen that, with high probability, the expected increment is on the
order of (cf. (32))

$\displaystyle\mathbb{E}\left[\langle\bm{u}^{1},\bm{u}^{\star}\rangle\mid\bm{u}^{0}\right]$
$\displaystyle-\langle\bm{u}^{0},\bm{u}^{\star}\rangle\approx\eta\langle\bm{u}^{0},\bm{u}^{\star}\rangle^{2}\|\bm{u}^{\star}\|_{2}^{2}\lesssim\frac{\eta\,\mathrm{poly}\log(d)}{d}\|\bm{u}^{\star}\|_{2}^{4}\|\bm{u}^{0}\|_{2}^{2},$
(33)

which could be quite small as it depends quadratically on the current
correlation $\langle\bm{u}^{0},\bm{u}^{\star}\rangle$.

If we were to hope that the favorable population-level analysis captures more
or less the finite-sample dynamics, we would need to ensure that the
variability of the gradient update is well-controlled. Towards this, let us
compute the variance of $\langle\bm{u}^{1},\bm{u}^{\star}\rangle$, assuming
that
$\frac{\|\bm{u}^{0}\|_{\infty}}{\|\bm{u}^{0}\|_{2}}\asymp\frac{\|\bm{u}^{\star}\|_{\infty}}{\|\bm{u}^{\star}\|_{2}}\asymp\frac{\mathrm{poly}\log(d)}{\sqrt{d}}$:

$\displaystyle\mathsf{Var}\left(\langle\bm{u}^{1},\bm{u}^{\star}\rangle\mid\bm{u}^{0}\right)$
$\displaystyle\asymp\mathsf{Var}\Bigg{(}\frac{\eta}{p}\sum_{1\leq j,k,l\leq
d}\big{(}\chi_{jkl}-p\big{)}\left(u_{j}^{0}u_{k}^{0}u_{l}^{0}-u_{j}^{\star}u_{k}^{\star}u_{l}^{\star}\right)u_{j}^{0}u_{k}^{0}u_{l}^{\star}\Bigg{)}$
$\displaystyle\asymp\frac{\eta^{2}}{p}\sum_{1\leq j,k,l\leq
d}\left(u_{j}^{0}u_{k}^{0}u_{l}^{0}+u_{j}^{\star}u_{k}^{\star}u_{l}^{\star}\right)^{2}\left(u_{j}^{0}u_{k}^{0}u_{l}^{\star}\right)^{2}\asymp\frac{\eta^{2}\mathrm{poly}\log(d)}{pd^{3}}\|\bm{u}^{\star}\|_{2}^{8}\|\bm{u}^{0}\|_{2}^{4}.$

In other words, the typical size of the variability of
$\langle\bm{u}^{1},\bm{u}^{\star}\rangle$ is about the order of
$\frac{\eta\mathrm{poly}\log(d)}{\sqrt{pd^{3}}}\|\bm{u}^{\star}\|_{2}^{4}\|\bm{u}^{0}\|_{2}^{2}$,
which dominates (in fact, is order-of-magnitudes larger than) the mean
increment (33) unless

$p\gtrsim\frac{\mathrm{poly}\log(d)}{d}.$ (34)

The sample size corresponding to (34) is, however, considerably larger than
the computation limit $p\asymp\frac{\mathrm{poly}\log(d)}{d^{1.5}}$. The
presence of a large variance implies highly volatile dynamics of randomly
initialized GD, thus casting doubt on its efficiency in the most challenging
sample-starved regime.

In summary, the main issue stems from the quadratic dependence of the expected
increment (33) on the correlation $\langle\bm{u}^{0},\bm{u}^{\star}\rangle$,
which can be exceedingly small if $\bm{u}^{0}$ is randomly initialized.

##### Initialization via the tensor power method (TPM).

Another alternative for initialization is the tensor power method, which has
recently gained popularity in the context of learning latent-variable models
[AGH+14, AGJ17]. Nevertheless, the TPM (with random initialization) suffers
from the same high-volatility issue as randomly initialized GD. The argument
for this would be nearly identical to the one presented above, and is hence
omitted. Instead, we invoke a perturbation analysis result in [AGH+14, Theorem
5.1] to illustrate the insufficiency of the TPM.

Recall that
$\frac{1}{p}\bm{T}=\bm{T}^{\star}+\big{(}\frac{1}{p}\bm{T}-\bm{T}^{\star}\big{)}$.
A critical issue is that the perturbation bound in [AGH+14, Theorem 5.1]
requires the tensor perturbation to be exceedingly small, namely,

$\big{\|}\tfrac{1}{p}\bm{T}-\bm{T}^{\star}\big{\|}\lesssim 1/d.$ (35)

This, however, cannot possibly hold if the sample size is merely
$p\asymp\frac{\mathrm{poly}\log(d)}{d^{1.5}}$ (in which case one only expects
a spectral norm bound on the order of
$\|p^{-1}\bm{T}-\bm{T}^{\star}\|\lesssim\frac{1}{\mathrm{poly}\log(d)}$ shown
in Corollary D.3 even in the absence of noise). In light of all this, existing
stability analysis of the TPM does not imply either sample efficiency or
computational efficiency.

##  4 Related work

One of the most natural ideas for solving tensor completion is to first unfold
the tensor data into matrices, followed by proper convex relaxation commonly
adopted for low-rank matrix completion. Given that there are more than one
ways to matricize a tensor, several prior work has explored the design of
matrix norms that can exploit the tensor structure more effectively [THK10,
GRY11, LMWY13, RPP13, LFC+16, MHWG14]. Such algorithms have been robustified
to enable reliable recovery against sparse outliers as well [GQ14]. For the
most part, however, such unfolding-based convex relaxation necessarily incur
loss of structural information, which is particularly severe when handling
odd-order tensors. The sample complexity developed for this paradigm is often
sub-optimal vis-a-vis the computational limits (namely, minimal sample
complexity achievable by polynomial-time algorithms).

Motivated by the above sub-optimality issue, [YZ16, YZ17] proposed to minimize
instead the tensor nuclear norm subject to data constraints, which provably
allows for reduced sample complexity. The issue, however, is that computing
the tensor nuclear norm itself is already computationally intractable, thus
limiting its applicability to even moderate-dimensional problems. Similar
findings have also been discovered for tensor atomic norm minimization
[DBBG19]. When restricted to polynomial-time algorithms, the best statistical
guarantees are often attained via convex relaxation tailored to the sum-of-
squares hierarchy [BM16]; the resulting computational cost, however, remains
prohibitively high for practical large-scale problems. Another matrix nuclear
norm minimization algorithm has been proposed based on promoting certain
structures on certain factor matrices [LSC+14]. Developing statistical
guarantees is, however, not the focal point of this work.

Moving beyond convex relaxation, a number of prior papers have developed
nonconvex algorithms for tensor completion, examples including iterative hard
thresholding [RSS17], alternating minimization [JO14, WAA16, XHYS15], tensor
SVD [ZA17], optimization on manifold [XY17, KM16, Ste16], proximal average
algorithm with nonconvex regularizer [Yao18], and block coordinate decent
[JHZ+16, XY13]. When it comes to the model considered herein, these algorithms
either lack optimal statistical guarantees, or come with a computational cost
that is significantly higher than a linear-time algorithm.

The algorithm and theory that we develop are largely inspired by the recent
advances of nonconvex optimization algorithms for low-rank matrix recovery
problems [KMO10a, KMO10b, CLS15, CC17, SL16, YPCC16, CW15]. The main
theoretical tool — the leave-one-out analysis — is a powerful technique that
has proved successful in various other statistical problems [EK15, CFMW19,
AFWZ17, MWCC17, ZB18, CCFM19, CCF+19, LZT19, CFMY19, DC18, PW19]. There are
several major differences between the analysis of nonconvex tensor completion
and that of nonconvex matrix recovery. For instance, our initialization scheme
is substantially more complicated than the matrix recovery counterpart, thus
requiring much more sophisticated analysis; in addition, the local convergence
stage of tensor completion does not suffer from rotational ambiguity (which
often appears in nonconvex matrix completion), and hence we only need to
handle permutational ambiguity.

In addition, the current paper focuses on non-adaptive uniform random
sampling. If there is freedom in designing the sampling mechanism, then one
can often expect improved performance; see [KS13, Zha19] as examples.
Fundamental criteria that enable perfect low-CP-rank tensor completion have
been studied in [AW17].

Tensor completion is simply a special example of the tensor recovery
literature. There is a large body of results tackling various other tensor
recovery and estimation problems, including but not limited to tensor
decomposition [Kol01, KB09, AGH+14, AGJ14, TS15, KOKC13, HSSS16, GHJY15,
ZKOM18, SDLF+17, SLLC17, GM17], tensor SVD and factorization [ZX18, KBHH13,
ZA17], and tensor regression and sketching [RSS17, HZC20, CRY19, HWW+19]. The
algorithmic ideas explored in this paper might have implications for these
tensor-related problems as well.

##  5 Analysis

In this section, we outline the proof of Theorem 2.8. The proof of Corollary
2.9 is deferred to Appendix C. The analysis is divided into three parts:

  * •

In Section 5.1, we show that given an initial estimate sufficiently close to
the ground truth, vanilla gradient descent converges linearly. These are
formalized in Lemmas 5.3 and 5.6.

  * •

Sections 5.2-5.3 provide statistical guarantees for the two steps of the
initialization procedure; see Theorems 5.9.

  * •

Under the assumptions of Theorem 2.8, one can see that the initialization
satisfies the requirement of linear convergence of vanilla gradient descent.
Therefore, Theorem 2.8 immediately follows from the results in Sections
5.1-5.3.

###  5.1 Analysis for local convergence of GD

In this section, we demonstrate that: if the initialization is reasonably
good, then vanilla gradient descent converges linearly to a solution with the
desired statistical accuracy. We postpone the analysis for initialization to
Sections 5.2-5.3 for convenience of presentation.

####  5.1.1 Preliminaries: gradient and Hessian calculation

First of all, using our notation $\times^{\mathsf{seq}}$ defined in (21), we
can write

$\displaystyle\nabla f\left(\bm{U}\right)$
$\displaystyle=\frac{1}{p}\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{1\leq
i\leq r}\bm{u}_{i}^{\otimes
3}-\bm{T}^{\star}-\bm{E}\Big{)}\times^{\mathsf{seq}}_{1}\bm{U}\times^{\mathsf{seq}}_{2}\bm{U}.$
(36)

Next, we find it convenient to define an auxiliary loss function
$f_{\mathsf{clean}}\left(\bm{U}\right):\mathbb{R}^{d\times
r}\to\mathbb{R}_{+}$ that corresponds to the noiseless case:

$\displaystyle
f_{\mathsf{clean}}\left(\bm{U}\right)=\frac{1}{6p}\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{1\leq
i\leq r}\bm{u}_{i}^{\otimes
3}-\bm{T}^{\star}\Big{)}\Big{\|}^{2}_{\mathrm{F}}.$ (37)

The gradient of $f_{\mathsf{clean}}$ w.r.t. $\bm{u}_{s}$ ($1\leq s\leq r$) is
thus given by

$\displaystyle\nabla_{\bm{u}_{s}}f_{\mathsf{clean}}\left(\bm{U}\right)$
$\displaystyle=\frac{1}{p}\,\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{1\leq
i\leq r}\bm{u}_{i}^{\otimes
3}-\bm{T}^{\star}\Big{)}\times_{1}\bm{u}_{s}\times_{2}\bm{u}_{s},\qquad 1\leq
s\leq r,$ (38)

and hence one can write

$\displaystyle\nabla f_{\mathsf{clean}}\left(\bm{U}\right)$
$\displaystyle=\frac{1}{p}\,\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{1\leq
i\leq r}\bm{u}_{i}^{\otimes
3}-\bm{T}^{\star}\Big{)}\times^{\mathsf{seq}}_{1}\bm{U}\times^{\mathsf{seq}}_{2}\bm{U}.$
(39)

This clearly satisfies

$\displaystyle\nabla f\left(\bm{U}\right)=\nabla
f_{\mathsf{clean}}\left(\bm{U}\right)-\frac{1}{p}\mathcal{P}_{\Omega}(\bm{E})\times^{\mathsf{seq}}_{1}\bm{U}\times^{\mathsf{seq}}_{2}\bm{U}.$
(40)

Moreover, direct algebraic manipulations give that: for any matrix
$\bm{V}=\left[\bm{v}_{1},\dots,\bm{v}_{r}\right]\in\mathbb{R}^{d\times r}$,

$\displaystyle\mathsf{vec}\left(\bm{V}\right)^{\top}\nabla^{2}f_{\mathsf{clean}}\left(\bm{U}\right)\mathsf{vec}\left(\bm{V}\right)$
$\displaystyle=\frac{1}{3p}\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{1\leq
s\leq
r}\bm{u}_{s}\otimes\bm{u}_{s}\otimes\bm{v}_{s}+\bm{u}_{s}\otimes\bm{v}_{s}\otimes\bm{u}_{s}+\bm{v}_{s}\otimes\bm{u}_{s}\otimes\bm{u}_{s}\Big{)}\Big{\|}_{\mathrm{F}}^{2}$
$\displaystyle\quad+\frac{2}{p}\Big{\langle}\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{s\in[r]}\bm{u}_{s}^{\otimes
3}-\bm{T}^{\star}\Big{)},\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{v}_{s}\otimes\bm{u}_{s}\Big{\rangle},$
(41)

where $\mathsf{vec}(\bm{V})$ denotes the vectorization of $\bm{V}$.

####  5.1.2 Local strong convexity and smoothness

At the heart of our analysis is a crucial geometric property of the objective
function, that is, the noiseless loss function $f_{\mathsf{clean}}$ behaves
like a locally strongly convex and smooth function. This fact, which is
formally stated in the following lemma, is the key enabler of fast local
convergence of vanilla GD.

######  Lemma 5.1 (Local strong convexity and smoothness).

Suppose that the sample complexity and the rank satisfy

$\displaystyle p\geq
c_{0}\max\left\\{\frac{\log^{3}d}{d^{3/2}},\frac{\mu^{2}r^{2}\log
d}{d^{2}}\right\\}\qquad r\leq c_{1}\sqrt{\frac{d}{\mu}}$ (42)

for some sufficiently large (resp. small) constant $c_{0}>0$ (resp.
$c_{1}>0$). Then with probability greater than $1-O(d^{-10})$,

$\displaystyle\frac{1}{2}\lambda^{\star
4/3}_{\min}\left\|\bm{V}\right\|_{\mathrm{F}}^{2}\leq\mathsf{vec}\left(\bm{V}\right)^{\top}\nabla^{2}f_{\mathsf{clean}}\left(\bm{U}\right)\,\mathsf{vec}\left(\bm{V}\right)\leq
4\lambda^{\star 4/3}_{\max}\left\|\bm{V}\right\|_{\mathrm{F}}^{2}$ (43)

holds simultaneously for all $\bm{V}\in\mathbb{R}^{d\times r}$ and all
$\bm{U}\in\mathbb{R}^{d\times r}$ obeying

$\|\bm{U}-\bm{U}^{\star}\|_{\mathrm{F}}\leq\delta\left\|\bm{U}^{\star}\right\|_{\mathrm{F}}\quad\text{and}\quad\left\|\bm{U}-\bm{U}^{\star}\right\|_{2,\infty}\leq\delta\left\|\bm{U}^{\star}\right\|_{2,\infty}.$
(44)

Here, $\delta\leq c_{2}/(\mu^{3/2}r)$ for some sufficiently small constant
$c_{2}>0$.

###### Proof.

See Appendix A.1. ∎

In order to invoke Lemma 5.1, one needs to make sure that the decision matrix
$\bm{U}$ of interest (e.g. $\bm{U}^{t}$ in the GD sequence) satisfies the
condition (44). This, however, is a fairly stringent condition, as it requires
$\bm{U}$ to be close to the truth in every single row.

####  5.1.3 Leave-one-out gradient descent sequences

Motivated by the analytical framework developed for low-rank matrix recovery
[MWCC17, CLL19], we introduce the following leave-one-out sequences, which
play a crucial role in guaranteeing that the entire trajectory
$\\{\bm{U}^{t}\\}_{t\geq 0}$ satisfies the condition (44) as required in Lemma
5.1.

Specifically, we define for each $1\leq m\leq d$ the following auxiliary loss
function:

$\displaystyle
f^{\left(m\right)}\left(\bm{U}\right)\triangleq\frac{1}{6p}\,\Big{\|}\mathcal{P}_{\Omega_{-m}}\Big{(}\sum\nolimits_{1\leq
s\leq r}\bm{u}_{s}^{\otimes
3}-\bm{T}^{\star}-\bm{E}\Big{)}\Big{\|}_{\mathrm{F}}^{2}+\frac{1}{6}\,\Big{\|}\mathcal{P}_{m}\Big{(}\sum\nolimits_{1\leq
s\leq r}\bm{u}_{s}^{\otimes
3}-\bm{T}^{\star}\Big{)}\Big{\|}_{\mathrm{F}}^{2},$ (45)

where

  * •

$\mathcal{P}_{\Omega_{m}}$: the projection onto the subspace of tensors
supported on $\\{(i,j,k)\in\Omega\colon i=m\text{ or }j=m\text{ or }k=m\\}$;

  * •

$\mathcal{P}_{\Omega_{-m}}$: the projection onto the subspace of tensors
supported on $\\{(i,j,k)\in\Omega\colon i\neq m\text{ and }j\neq m\text{ and
}k\neq m\\}$;

  * •

$\mathcal{P}_{m}$: the projection onto the subspace of tensors supported on
$\\{(i,j,k)\in[d]^{3}\colon i=m\text{ or }j=m\text{ or }k=m\\}$.

In words, this function is obtained by replacing all data at locations
$\\{(i,j,k)\in[d]^{3}\colon i=m\text{ or }j=m\text{ or }k=m\\}$ by their
expected values, thus removing all randomness associated with this location
subset. The gradient of $f^{\left(m\right)}(\bm{U})$ w.r.t. $\bm{u}_{s}$
($1\leq s\leq r$) can be computed as:

$\displaystyle\begin{split}\nabla_{\bm{u}_{s}}f^{\left(m\right)}(\bm{U})&=\frac{1}{p}\mathcal{P}_{\Omega_{-m}}\Big{(}\sum\nolimits_{1\leq
s\leq r}\bm{u}_{s}^{\otimes
3}-\bm{T}^{\star}-\bm{E}\Big{)}\times_{1}\bm{u}_{s}\times_{2}\bm{u}_{s}+\mathcal{P}_{m}\Big{(}\sum\nolimits_{1\leq
s\leq r}\bm{u}_{s}^{\otimes
3}-\bm{T}^{\star}\Big{)}\times_{1}\bm{u}_{s}\times_{2}\bm{u}_{s}.\end{split}$
(46)

We then denote by $\big{\\{}\bm{U}^{t,(m)}\big{\\}}_{t\geq 0}$ the iterative
sequence obtained by running gradient descent w.r.t. the leave-one-out loss
$f^{\left(m\right)}(\cdot)$; see Algorithm 4. By construction, as long as
$\bm{U}^{0,\left(m\right)}$ is independent of the sampling locations and the
noise associated with the locations $\\{(i,j,k)\in\Omega\colon i=m\text{ or
}j=m\text{ or }k=m\\}$ (which holds true as detailed momentarily), then the
entire trajectory $\big{\\{}\bm{U}^{t,(m)}\big{\\}}_{t\geq 0}$ becomes
statistically independent of such randomness. This is a crucial property that
allows us to decouple the complicated statistical dependency.

Algorithm 4 The $m$-th leave-one-out sequence

1:Generate an initial estimate $\bm{U}^{0,\left(m\right)}$ via Algorithm 5.

2:for $t=0,1,\dots,t_{0}-1$ do

3: $\bm{U}^{t+1,\left(m\right)}=\bm{U}^{t,\left(m\right)}-\eta_{t}\nabla
f^{\left(m\right)}\big{(}\bm{U}^{t,\left(m\right)}\big{)}$.

####  5.1.4 Key lemmas

The proof for local linear convergence of GD is inductive in nature, which
proceeds on the basis of the following set of inductive hypotheses. As we
shall see in Corollary 5.11 in Section 5.3, this set of inductive hypotheses —
modulo some global permutation — is valid with high probability when $t=0$. In
order to simplify presentation, we remove the consideration of the global
permutation factor throughout this section (namely, we assume that the
following holds for $\bm{U}^{0}\bm{\Pi}^{0}$ with some permutation matrix
$\bm{\Pi}^{0}\in\mathbb{R}^{r\times r}$ obeying $\bm{\Pi}^{0}=\bm{I}$. Our
inductive hypotheses are summarized as follows:

Key hypotheses for the gradient update stage:

$\displaystyle\big{\|}\bm{U}^{t}-\bm{U}^{\star}\big{\|}_{\mathrm{F}}$
$\displaystyle\leq\left(C_{1}\rho^{t}\mathcal{E}_{\mathsf{local}}+C_{2}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{\mathrm{F}};$ (47a)
$\displaystyle\big{\|}\bm{U}^{t}-\bm{U}^{\star}\big{\|}_{2,\infty}$
$\displaystyle\leq\left(C_{3}\rho^{t}\mathcal{E}_{\mathsf{local}}+C_{4}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{2,\infty};$ (47b)
$\displaystyle\big{\|}\bm{U}^{t}-\bm{U}^{t,\left(m\right)}\big{\|}_{\mathrm{F}}$
$\displaystyle\leq\left(C_{5}\rho^{t}\mathcal{E}_{\mathsf{local}}+C_{6}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{2,\infty};$ (47c)
$\displaystyle\left\|\big{(}\bm{U}^{t,\left(m\right)}-\bm{U}^{\star}\big{)}_{m,:}\right\|_{2}$
$\displaystyle\leq\left(C_{7}\rho^{t}\mathcal{E}_{\mathsf{local}}+C_{8}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{2,\infty};$ (47d)

for some quantity $\mathcal{E}_{\mathsf{local}}>0$ (depending possibly on
$\mu$ and $r$) and some constants $C_{1},\cdots,C_{8}>0$. These exist a few
straightforward consequences of the hypotheses (47), which we record in the
following lemma.

######  Lemma 5.2.

Assume that the hypotheses (47) hold, then we have

$\displaystyle\big{\|}\bm{U}^{t,\left(m\right)}-\bm{U}^{\star}\big{\|}_{\mathrm{F}}$
$\displaystyle\leq\left(2C_{1}\rho^{t}\mathcal{E}_{\mathsf{local}}+2C_{2}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{\mathrm{F}},$ (48)
$\displaystyle\big{\|}\bm{U}^{t,\left(m\right)}-\bm{U}^{\star}\big{\|}_{2,\infty}$
$\displaystyle\leq\left(\left(C_{3}+C_{5}\right)\rho^{t}\mathcal{E}_{\mathsf{local}}+\left(C_{4}+C_{6}\right)\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{2,\infty}.$ (49)

###### Proof.

See Appendix A.2. ∎

Our proof for the hypotheses (47) is inductive in nature: we would like to
show that if the hypotheses in (47) hold for the $t$-th iteration, then they
continue to be valid for the $(t+1)$-th iteration. We shall justify each of
the above hypotheses inductively through the following lemmas.

######  Lemma 5.3.

Suppose that

$\displaystyle p\geq
c_{0}\frac{\mu^{3}r^{2}\log^{3}d}{d^{3/2}},\quad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{\mu^{3/2}r\sqrt{d\log d}},\quad\text{and}\quad r\leq
c_{2}\sqrt{\frac{d}{\mu}}$

for some sufficiently large constant $c_{0}>0$ and some sufficiently small
constant $c_{1},c_{2}>0$. Assume that the hypotheses (47) hold for the $t$-th
iteration and $\mathcal{E}_{\mathsf{local}}\leq
c_{3}/\big{(}\mu^{3/2}r\big{)}$ for some sufficiently small constant
$c_{3}>0$. Then with probability at least $1-O(d^{-10})$,

$\displaystyle\left\|\bm{U}^{t+1}-\bm{U}^{\star}\right\|_{\mathrm{F}}\leq\left(C_{1}\rho^{t+1}\mathcal{E}_{\mathsf{local}}+C_{2}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{\mathrm{F}},$ (50)

provided that $0<\eta\leq\lambda_{\min}^{\star
4/3}/\big{(}32\lambda_{\max}^{\star 8/3}\big{)}$,
$1-\big{(}\lambda_{\min}^{\star 4/3}/5\big{)}\eta\leq\rho<1$, and $C_{2}$ is
sufficiently large.

###### Proof.

See Appendix A.3. ∎

######  Lemma 5.4.

Suppose that

$\displaystyle p\geq
c_{0}\frac{\mu^{3}r^{2}\log^{3}d}{d^{3/2}},\quad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{\mu^{3/2}r\sqrt{d\log d}},\quad\text{and}\quad r\leq
c_{2}\sqrt{\frac{d}{\mu}}$

for some sufficiently large constant $c_{0}>0$ and some sufficiently small
constant $c_{1},c_{2}>0$. Assume that the hypotheses (47) hold for the $t$-th
iteration and $\mathcal{E}_{\mathsf{local}}\leq
c_{3}/\big{(}\mu^{3/2}r\big{)}$ for some sufficiently small constant
$c_{3}>0$. Then with probability at least $1-O(d^{-10})$, one has

$\displaystyle\big{\|}\bm{U}^{t+1,\left(m\right)}-\bm{U}^{t+1}\big{\|}_{\mathrm{F}}\leq\left(C_{5}\rho^{t+1}\mathcal{E}_{\mathsf{local}}+C_{6}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{2,\infty},$ (51)

provided that $0<\eta\leq\lambda_{\min}^{\star
4/3}/\big{(}32\lambda_{\max}^{\star 8/3}\big{)}$,
$1-\big{(}\lambda_{\min}^{\star 4/3}/5\big{)}\eta\leq\rho<1$ and $C_{6}$ is
sufficiently large.

###### Proof.

See Appendix A.4. ∎

######  Lemma 5.5.

Suppose that

$\displaystyle p\geq
c_{0}\frac{\mu^{3}r^{2}\log^{3}d}{d^{3/2}},\quad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{\mu^{3/2}r\sqrt{d\log d}},\quad\text{and}\quad r\leq
c_{2}\sqrt{\frac{d}{\mu}}$

for some sufficiently large constant $c_{0}>0$ and some sufficiently small
constant $c_{1},c_{2}>0$. Assume that the hypotheses (47) hold for the $t$-th
iteration and $\mathcal{E}_{\mathsf{local}}\leq
c_{3}/\big{(}\mu^{3/2}r\big{)}$ for some sufficiently small constant
$c_{3}>0$. Then with probability at least $1-O(d^{-10})$, one has

$\displaystyle\left\|\big{(}\bm{U}^{t+1,\left(m\right)}-\bm{U}^{\star}\big{)}_{m,:}\right\|_{2}\leq\left(C_{7}\rho^{t+1}\mathcal{E}_{\mathsf{local}}+C_{8}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{2,\infty},$ (52)

provided that $0<\eta\leq\lambda_{\min}^{\star
4/3}/\big{(}32\lambda_{\max}^{\star 8/3}\big{)}$,
$1-\big{(}\lambda_{\min}^{\star 4/3}/5\big{)}\eta\leq\rho<1$, $C_{7}$ and
$C_{8}$ are sufficiently large.

###### Proof.

See Appendix A.5. ∎

######  Lemma 5.6.

Suppose that

$\displaystyle p\geq
c_{0}\frac{\mu^{3}r^{2}\log^{3}d}{d^{3/2}},\quad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{\mu^{3/2}r\sqrt{d\log d}},\quad\text{and}\quad r\leq
c_{2}\sqrt{\frac{d}{\mu}}$

for some sufficiently large constant $c_{0}>0$ and some sufficiently small
constant $c_{1},c_{2}>0$. Assume that the hypotheses (47) hold for the $t$-th
iteration and $\mathcal{E}_{\mathsf{local}}\leq
c_{3}/\big{(}\mu^{3/2}r\big{)}$ for some sufficiently small constant
$c_{3}>0$. Then with probability at least $1-O(d^{-10})$, one has

$\displaystyle\left\|\bm{U}^{t+1}-\bm{U}^{\star}\right\|_{2,\infty}\leq\left(C_{3}\rho^{t+1}\mathcal{E}_{\mathsf{local}}+C_{4}\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log
d}{p}}\right)\left\|\bm{U}^{\star}\right\|_{2,\infty},$ (53)

provided that $0<\eta\leq\lambda_{\min}^{\star
4/3}/\big{(}32\lambda_{\max}^{\star 8/3}\big{)}$,
$1-\big{(}\lambda_{\min}^{\star 4/3}/5\big{)}\eta\leq\rho<1$,
$C_{3}/\left(C_{5}+C_{7}\right)$ and $C_{4}/\left(C_{6}+C_{8}\right)$ are both
sufficiently large.

###### Proof.

See Appendix A.6. ∎

The proofs of the above key lemmas are postponed to Appendix A.

###  5.2 Analysis for initialization: Part 1 (subspace estimation)

####  5.2.1 Key results

The aim of this subsection is to demonstrate that the subspace estimate
$\bm{U}$ computed by Algorithm 2 is sufficiently close to the space spanned by
the true tensor factors. Given that the columns of
$\bm{U}^{\star}=[\bm{u}_{1}^{\star},\cdots,\bm{u}_{r}^{\star}]$ are in general
not orthogonal to each other, we shall define
$\bm{U}_{\mathsf{orth}}^{\star}\in\mathbb{R}^{d\times r}$ as follows (obtained
by proper orthonormalization) :

$\displaystyle\bm{U}_{\mathsf{orth}}^{\star}:={\bm{U}}^{\star}\big{(}{\bm{U}}^{\star\top}{\bm{U}}^{\star}\big{)}^{-\frac{1}{2}}.\qquad$
(54)

This matrix $\bm{U}_{\mathsf{orth}}^{\star}$ reflects the rank-$r$ principal
subspace of
$\bm{A}^{\star}\bm{A}^{\star\top}=\sum_{i}\|\bm{u}_{i}^{\star}\|_{2}^{4}\bm{u}_{i}^{\star}\bm{u}_{i}^{\star\top}$,
where we recall that $\bm{A}^{\star}\in\mathbb{R}^{d\times d^{2}}$ is the
mode-1 matricization of $\bm{T}^{\star}$. In addition, we define the rotation
matrix

$\displaystyle\bm{R}:=\operatorname*{arg\,min}_{\bm{Q}\in\mathcal{O}^{r\times
r}}\leavevmode\nobreak\
\left\|\bm{U}\bm{Q}-\bm{U}^{\star}_{\mathsf{orth}}\right\|_{\mathrm{F}},$ (55)

where $\mathcal{O}^{r\times r}$ stands for the set of $r\times r$ orthonormal
matrices. This can be viewed as the global rotation matrix that best aligns
the two subspaces represented by $\bm{U}$ and $\bm{U}^{\star}_{\mathsf{orth}}$
respectively.

Equipped with the above notation, we can invoke [CLC+20, Corollary 1] to
arrive at the following lemma, which upper bounds the distance between our
subspace estimate $\bm{U}$ and the ground truth
$\bm{U}^{\star}_{\mathsf{orth}}$.

######  Lemma 5.7.

There exist some universal constants $c_{0},c_{1},c_{2}>0$ such that if

$\displaystyle p\geq
c_{0}\frac{\mu^{2}r\log^{2}d}{d^{3/2}},\quad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{d^{3/4}\sqrt{\log d}},\quad\text{and}\quad r\leq
c_{2}\sqrt{\frac{d}{\mu}},$

then with probability $1-O\left(d^{-10}\right)$, the subspace estimate
$\bm{U}$ computed by Algorithm 2 obeys

$\displaystyle\left\|\bm{U}\bm{R}-\bm{U}_{\mathsf{orth}}^{\star}\right\|$
$\displaystyle\lesssim\mathcal{E}_{\mathsf{se}},$ (56a)
$\displaystyle\left\|\bm{U}\bm{R}-\bm{U}_{\mathsf{orth}}^{\star}\right\|_{2,\infty}$
$\displaystyle\lesssim\mathcal{E}_{\mathsf{se}}\sqrt{\frac{\mu r}{d}},$ (56b)

where $\bm{U}_{\mathsf{orth}}^{\star}$ and $\bm{R}$ are defined respectively
in (54) and (55), and

$\displaystyle\mathcal{E}_{\mathsf{se}}:=\frac{\mu^{2}r\log
d}{d^{3/2}p}+\sqrt{\frac{\mu^{2}r\log
d}{d^{2}p}}+\frac{\sigma^{2}}{\lambda_{\min}^{\star 2}}\frac{d^{3/2}\log
d}{p}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log d}{p}}+\frac{\mu
r}{d}.$ (57)

In a nutshell, Lemma 5.7 asserts that: under our sample size, noise and rank
conditions, Algorithm 2 produces reliable estimates of the subspace spanned by
the low-rank tensor factors $\\{\bm{u}_{i}^{\star}\\}_{1\leq i\leq r}$. The
theorem quantifies the subspace distance in terms of both the spectral norm
and $\|\cdot\|_{2,\infty}$, where the latter bound often reflects a
considerably stronger sense of proximity compared to the former one.

As it turns out, in order to facilitate analysis for the subsequent stages, we
need to introduce certain leave-one-out sequences as well, which we detail in
the next subsection.

####  5.2.2 Leave-one-out sequences for subspace estimation

The key idea of the leave-one-out analysis is to create auxiliary leave-one-
out sequences that are (1) independent of a small fraction of the data; (2)
sufficiently close to the true estimates. We introduce the following auxiliary
tensor and $d\times d^{2}$-dimensional matrix for each $1\leq m\leq d$:

$\displaystyle\bm{T}^{\left(m\right)}$
$\displaystyle:=\mathcal{P}_{\Omega_{-m}}\left(\bm{T}\right)+p\mathcal{P}_{m}\left(\bm{T}^{\star}\right)\in\mathbb{R}^{d\times
d\times d},$ (58) $\displaystyle\bm{A}^{\left(m\right)}$
$\displaystyle:=\text{mode-1 matricization of
}\tfrac{1}{p}\bm{T}^{\left(m\right)}.$ (59)

By construction, $\bm{T}^{\left(m\right)}$ and $\bm{A}^{\left(m\right)}$ are
independent of $\mathcal{P}_{\Omega_{m}}\left(\bm{E}\right)$, where we recall
that

$\displaystyle\Omega_{-m}$ $\displaystyle:=\\{(i,j,k)\in\Omega\colon i\neq
m\text{ and }j\neq m\text{ and }k\neq m\\},$ (60) $\displaystyle\Omega_{m}$
$\displaystyle:=\\{(i,j,k)\in\Omega\colon i=m\text{ or }j=m\text{ or }k=m\\}.$
(61)

We are now ready to introduce the auxiliary leave-one-out procedure for
subspace estimation. Similar to the matrix $\bm{B}$ in Algorithm 2 (whose
eigenspace serves as an estimate of the column space of $\bm{U}^{\star}$), we
define an auxiliary matrix $\bm{B}^{\left(m\right)}\in\mathbb{R}^{d\times d}$
as follows:

$\displaystyle\bm{B}^{\left(m\right)}=\mathcal{P}_{\mathsf{off}\text{-}\mathsf{diag}}\big{(}\bm{A}^{\left(m\right)}\bm{A}^{\left(m\right)\top}\big{)},$
(62)

where $\mathcal{P}_{\mathsf{off}\text{-}\mathsf{diag}}(\cdot)$ (as already
defined in Section 3.1) extracts out off-diagonal entries from a matrix. The
rationale is simple: it can be easily verified that

$\displaystyle\mathbb{E}\big{[}\bm{B}^{\left(m\right)}\big{]}=\bm{B}^{\star}-\mathcal{P}_{\mathsf{diag}}\left(\bm{B}^{\star}\right),\qquad\bm{B}^{\star}:=\bm{A}^{\star}\bm{A}^{\star\top},$
(63)

where $\mathcal{P}_{\mathsf{diag}}(\cdot)$ extracts out the diagonal entries
of the matrix. This gives hope that the eigenspace of
$\bm{B}^{\left(m\right)}$ is also a reliable estimate of the column space of
$\bm{U}^{\star}$, provided that the diagonal entries of $\bm{B}^{\star}$ are
sufficiently small. Consequently, we shall compute
$\bm{U}^{0,\left(m\right)}\in\mathbb{R}^{d\times r}$ — a matrix whose columns
are the top-$r$ leading eigenvectors of $\bm{B}^{\left(m\right)}$. The
procedure is summarized in Algorithm 5.

Algorithm 5 The $m$-th leave-one-out sequence for spectral initialization

1:Let
$\bm{U}^{\left(m\right)}\bm{\Lambda}^{\left(m\right)}\bm{U}^{\left(m\right)\top}$
be the rank-$r$ eigen-decomposition of $\bm{B}^{\left(m\right)}$ defined in
(62).

2:Generate the initial estimate
$\bm{U}^{0,\left(m\right)}\in\mathbb{R}^{d\times r}$ from
$\bm{U}^{\left(m\right)}\in\mathbb{R}^{d\times r}$ using Algorithm 6.

The following lemma plays a crucial role in our analysis, which formalizes the
fact that the leave-one-out version $\bm{U}^{\left(m\right)}$ obtained by
Algorithm 5 is extremely close to $\bm{U}$.

######  Lemma 5.8.

There exist some universal constants $c_{0},c_{1},c_{2}>0$ such that if

$\displaystyle p\geq
c_{0}\frac{\mu^{2}r\log^{2}d}{d^{3/2}},\quad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{d^{3/4}\sqrt{\log d}},\quad\text{and}\quad r\leq
c_{2}\sqrt{\frac{d}{\mu}},$

then with probability $1-O\left(d^{-10}\right)$, the subspace estimate
$\bm{U}^{(m)}$ computed by Algorithm 5 obeys

$\displaystyle\big{\|}\bm{U}\bm{U}^{\top}-\bm{U}^{\left(m\right)}\bm{U}^{\left(m\right)\top}\big{\|}_{\mathrm{F}}$
$\displaystyle\lesssim\mathcal{E}_{\mathsf{loo}}\sqrt{\frac{\mu r}{d}}$ (64)

simultaneously for all $1\leq m\leq d$, where

$\displaystyle\mathcal{E}_{\mathsf{loo}}:=\frac{\mu^{2}r\log
d}{d^{3/2}p}+\sqrt{\frac{\mu^{2}r\log
d}{d^{2}p}}+\frac{\sigma^{2}}{\lambda_{\min}^{\star 2}}\frac{d^{3/2}\log
d}{p}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log d}{p}}.$ (65)

Lemma 5.8 follows immediately from the analysis of [CLC+20, Lemma 4]. As a
remark, the construction of the leave-one-out sequences herein is slightly
different from the one in [CLC+20]. However, it is straightforward to adapt
the proof of [CLC+20] to the case considered herein. We therefore omit the
proof for the sake of brevity.

###  5.3 Analysis for initialization: Part 2 (retrieval of individual tensor
factors)

####  5.3.1 Main results and leave-one-out sequences

This section justifies that the procedure presented in Algorithm 3 allows to
disentangle the tensor factors. For notational simplicity, we let

$\displaystyle\overline{\bm{u}}_{i}^{\star}:=\bm{u}_{i}^{\star}\,/\left\|\bm{u}_{i}^{\star}\right\|_{2},\qquad\lambda_{i}^{\star}:=\left\|\bm{u}_{i}^{\star}\right\|_{2}^{3},\qquad
1\leq i\leq d.$ (66)

Our result is this:

######  Theorem 5.9.

Fix any arbitrary small constant $\delta>0$. Assume that

$\displaystyle p\geq
c_{0}\frac{\mu^{2}r^{4}\log^{4}d}{d^{3/2}},\quad\frac{\sigma}{\lambda_{\min}^{\star}}\leq
c_{1}\frac{\sqrt{p}}{r^{3/2}d^{3/4}\log^{2}d},\quad r\leq
c_{2}\left(\frac{d}{\mu^{6}\log^{6}d}\right)^{1/6},$ $\displaystyle
L=c_{3}r^{2\kappa^{2}}\log^{3/2}r,\quad\epsilon_{\mathsf{th}}=c_{4}\Bigg{\\{}\frac{\mu
r\log
d}{d\sqrt{p}}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{rd\log^{2}d}{p}}+\sqrt{\frac{\mu
r\log d}{d}}\Bigg{\\}}$ (67)

for some sufficiently large universal constant $c_{0},c_{3}>0$ and some
sufficiently small universal constants $c_{1},c_{2},c_{4}>0$. Then with
probability exceeding $1-\delta$, there exists a permutation
$\pi(\cdot):[d]\mapsto[d]$ such that for all $1\leq i\leq r$, the tensor
factors $\\{\bm{w}^{i}\\}_{i=1}^{r}$ returned by Algorithm 3 satisfy

$\displaystyle\big{\|}\bm{w}^{i}-\overline{\bm{u}}^{\star}_{\pi(i)}\big{\|}_{2}$
$\displaystyle\lesssim\frac{\mu r\log
d}{d\sqrt{p}}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{rd\log^{2}d}{p}}+\sqrt{\frac{\mu
r\log d}{d}};$ (68a)
$\displaystyle\big{\|}\bm{w}^{i}-\overline{\bm{u}}^{\star}_{\pi(i)}\big{\|}_{\infty}$
$\displaystyle\lesssim\Bigg{\\{}\frac{\mu^{2}r\log^{4}d}{d^{3/2}p}+\frac{\mu
r\log^{3}d}{d\sqrt{p}}+\frac{\sigma^{2}}{\lambda_{\min}^{\star
2}}\frac{d^{3/2}\log^{4}d}{p}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{rd\log^{6}d}{p}}+\sqrt{\frac{\mu
r\log^{2}d}{d}}\Bigg{\\}}\sqrt{\frac{\mu r}{d}};$ (68b)
$\displaystyle\big{|}\lambda_{i}-\lambda_{\pi(i)}^{\star}\big{|}$
$\displaystyle\lesssim\Bigg{\\{}\frac{\mu r\log
d}{d\sqrt{p}}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{rd\log^{2}d}{p}}+\sqrt{\frac{\mu
r\log d}{d}}\Bigg{\\}}\lambda^{\star}_{\pi(i)}.$ (68c)

In short, this theorem asserts that the estimates returned by Algorithm 3 are
— up to global permutation — reasonably close to the ground truth under our
sample size and noise conditions. In order to establish this theorem and in
order to provide initial guesses for the leave-one-out GD sequences, we need
to produce a leave-one-out sequence tailored to this part of the algorithm.
Such auxiliary sequences are generated in a similar spirit as the previous
ones, and we summarize them in Algorithm 6. As usual, the resulting leave-one-
out estimates
$\big{\\{}\lambda^{\left(m\right)}_{i},\bm{w}^{i,\left(m\right)}\big{\\}}_{i=1}^{r}$
are statistically independent of
$\mathcal{P}_{\Omega_{m}}\left(\bm{E}\right)$.

In what follows, we gather a few key properties of the leave-one-out
estimates, which play a crucial role in the analysis.

Algorithm 6 The $m$-th leave-one-out sequence for retrieving individual tensor
components

1:Input: restart number $L$, threshold $\epsilon_{\mathsf{th}}$, subspace
estimate $\bm{U}^{\left(m\right)}\in\mathbb{R}^{d\times r}$ given by Algorithm
5.

2:for $\tau=1,\dots,L$ do

3: Recall the Gaussian vector $\bm{g}^{\tau}\sim\mathcal{N}(0,\bm{I}_{d})$
generated in Algorithm 3.

4:
$\big{(}\bm{\nu}^{\tau,\left(m\right)},\lambda_{\tau}^{\left(m\right)},\mathsf{spec}\text{-}\mathsf{gap}_{\tau}^{\left(m\right)}\big{)}\leftarrow\textsc{Retrieve-
one-tensor-
factor}(\bm{T}^{\left(m\right)},p,\bm{U}^{\left(m\right)},\bm{g}^{\tau})$.

5:Generate tensor factor estimates
$\displaystyle\big{\\{}\big{(}\bm{w}^{1,\left(m\right)},\lambda_{1}^{\left(m\right)}),\dots,(\bm{w}^{r,\left(m\right)},\lambda_{r}^{\left(m\right)}\big{)}\big{\\}}\leftarrow\textsc{Prune}(\big{\\{}\big{(}\bm{\nu}^{\tau,\left(m\right)},\lambda_{\tau}^{\left(m\right)},\mathsf{spec}\text{-}\mathsf{gap}_{\tau}^{\left(m\right)}\big{)}\big{\\}}_{\tau=1}^{L},\epsilon_{\mathsf{th}}).$

6:Output: an initial estimate
$\bm{U}^{0,\left(m\right)}=\big{[}\big{(}\lambda_{1}^{\left(m\right)}\big{)}^{1/3}\bm{w}^{1,\left(m\right)},\dots,\big{(}\lambda_{r}^{\left(m\right)}\big{)}^{1/3}\bm{w}^{r,\left(m\right)}\big{]}.$

######  Theorem 5.10.

Fix any arbitrarily small constant $\delta>0$. Instate the assumptions in
Theorem 5.9. With probability exceeding $1-\delta$, the permutation function
stated in Theorem 5.9 obeys that: for all $1\leq i\leq r$ and all $1\leq m\leq
d$:

$\displaystyle\big{\|}\bm{w}^{i}-\bm{w}^{i,\left(m\right)}\big{\|}_{2}$
$\displaystyle\lesssim\Bigg{\\{}\frac{\mu^{2}r\log^{3/2}d}{d^{3/2}p}+\frac{\mu\sqrt{r}\,\log
d}{d\sqrt{p}}+\frac{\sigma^{2}}{\lambda_{\min}^{\star
2}}\frac{d^{3/2}\log^{3/2}d}{p}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log^{2}d}{p}}\Bigg{\\}}\sqrt{\frac{\mu
r}{d}};$ (69a)
$\displaystyle\big{|}\lambda_{i}-\lambda_{i}^{\left(m\right)}\big{|}$
$\displaystyle\lesssim\Bigg{\\{}\frac{\mu^{2}r\log^{3/2}d}{d^{3/2}p}+\frac{\mu\sqrt{r}\,\log
d}{d\sqrt{p}}+\frac{\sigma^{2}}{\lambda_{\min}^{\star
2}}\frac{d^{3/2}\log^{3/2}d}{p}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log^{2}d}{p}}\Bigg{\\}}\sqrt{\frac{\mu
r}{d}}\,\lambda_{\max}^{\star};$ (69b)
$\displaystyle\big{|}\big{(}\bm{w}^{i,\left(m\right)}-\overline{\bm{u}}^{\star}_{\pi(i)}\big{)}_{m}\big{|}$
$\displaystyle\lesssim\Bigg{\\{}\frac{\sqrt{\mu
r}\log^{7/2}d}{d^{3/2}p}+\frac{\mu
r\log^{3}d}{d\sqrt{p}}+\frac{\sigma}{\lambda_{\min}^{\star}}\frac{\log^{4}d}{p}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{rd\log^{6}d}{p}}+\sqrt{\frac{\mu
r\log^{2}d}{d}}\Bigg{\\}}\sqrt{\frac{\mu r}{d}}.$ (69c)

With Theorems 5.9-5.10 in place, we can immediately establish a few desired
properties (particularly those specified in Section 5.1) of our initial
estimate, as asserted in the following corollary.

######  Corollary 5.11.

Fix any arbitrarily small constant $\delta>0$. Instate the assumptions in
Theorem 2.8. With probability exceeding $1-\delta$, the estimates $\bm{U}^{0}$
and $\bm{U}^{0,\left(m\right)}$ returned by Algorithm 3 and Algorithm 6
respectively satisfy the hypotheses (47) for $t=0$.

###### Proof.

See Appendix B.10. ∎

####  5.3.2 Analysis

Before we start with the proof, we first state the main idea. For the sake of
clarify, we define

$\displaystyle\bm{\theta}^{\tau}$
$\displaystyle:=\bm{U}\bm{U}^{\top}\bm{g}^{\tau},$ (70a)
$\displaystyle\bm{\theta}^{\tau,(m)}$
$\displaystyle:=\bm{U}^{(m)}\bm{U}^{(m)\top}\bm{g}^{\tau},$ (70b)
$\displaystyle\bm{M}^{\tau}$
$\displaystyle:=p^{-1}\bm{T}\times_{3}\bm{\theta}^{\tau},$ (70c)
$\displaystyle\bm{M}^{\tau,(m)}$
$\displaystyle:=p^{-1}\bm{T}^{(m)}\times_{3}\bm{\theta}^{\tau,(m)}.$ (70d)

In addition, let $\bm{\nu}^{\tau}$ be the top singular vector of
$\bm{M}^{\tau}$ obeying $\langle\bm{T},(\bm{\nu}^{\tau})^{\otimes
3}\rangle\geq 0$, and $\bm{\nu}^{\tau,(m)}$ the top singular vector of
$\bm{M}^{\tau,(m)}$ obeying
$\langle\bm{T}^{(m)},(\bm{\nu}^{\tau,(m)})^{\otimes 3}\rangle\geq 0$. Set

$\displaystyle\lambda_{\tau}:=\langle p^{-1}\bm{T},(\bm{\nu}^{\tau})^{\otimes
3}\rangle\qquad\text{and}\qquad\lambda_{\tau}^{(m)}:=\langle
p^{-1}\bm{T}^{(m)},(\bm{\nu}^{\tau,(m)})^{\otimes 3}\rangle.$ (71)

These are all computed in the function $\textsc{Retrieve-one-tensor-factor}($)
in the $\tau$-th round.

  1. 1.

We first show that for each $1\leq i\leq r$, there exists at least one trial
$1\leq\tau\leq L$ such that the $i$-th tensor factor
$\overline{\bm{u}}_{i}^{\star}$ is the top singular vector of the population
version of $\bm{T}\times_{3}\bm{\theta}^{\tau}$ (with respect to the missing
data and noise). In addition, the spectral gap is large enough to guarantee
accurate estimates.

  2. 2.

Next, we prove that given this spectral gap, the top singular vector
$\bm{\nu}^{\tau}$ of $\bm{T}\times_{3}\bm{\theta}^{\tau}$ is close to
$\overline{\bm{u}}^{\star}_{i}$ both in the $\|\cdot\|_{2}$ and
$\|\cdot\|_{\infty}$ norm. This also enables us to accurately estimate the
magnitude of $\bm{u}^{\star}_{i}$.

  3. 3.

Finally, we need to show that one can find those reliable estimates among $L$
random restarts. Combining the spectral gap information with the incoherence
condition that tensor components are nearly orthogonal to each other, our
selection procedure is guaranteed to recover all tensor factors.

Now we proceed to the proof. Without loss of generality, we prove the case for
$i=1$ in the sequel, i.e. there exists some $\tau\in\left[L\right]$ such that
$\bm{\nu}^{\tau}$ accurately recovers $\bm{u}^{\star}_{1}$. Together with the
union bound, this shows that we can find reliable estimates for all tensor
factors. We then conclude the proof by showing that Algorithm 3 is able to
find all of them without duplicates.

To this end, we find it convenient to introduce an auxiliary vector
$\bm{\gamma}^{\star\tau}=\left[\gamma_{1}^{\star\tau},\cdots,\gamma_{r}^{\star\tau}\right]^{\top}\in\mathbb{R}^{r}$
and its leave-one-out versions
$\bm{\gamma}^{\star\tau,\left(m\right)}=\big{[}\gamma_{1}^{\star\tau,\left(m\right)},\cdots,\gamma_{r}^{\star\tau,\left(m\right)}\big{]}^{\top}$
($1\leq m\leq d$) for each $1\leq\tau\leq L$ as follows:

$\displaystyle\gamma_{i}^{\star\tau}$
$\displaystyle:=\|\bm{u}_{i}^{\star}\|_{2}^{2}\,\big{\langle}\bm{u}_{i}^{\star},\bm{\theta}^{\tau}\big{\rangle}=\lambda_{i}^{\star}\big{\langle}\overline{\bm{u}}_{i}^{\star},\bm{\theta}^{\tau}\big{\rangle},$
(72a) $\displaystyle\gamma_{i}^{\star\tau,(m)}$
$\displaystyle:=\|\bm{u}_{i}^{\star}\|_{2}^{2}\,\big{\langle}\bm{u}_{i}^{\star},\bm{\theta}^{\tau,(m)}\big{\rangle}=\lambda_{i}^{\star}\big{\langle}\overline{\bm{u}}_{i}^{\star},\bm{\theta}^{\tau,(m)}\big{\rangle},$
(72b)

where $\overline{\bm{u}}_{i}^{\star}$ and $\lambda_{i}^{\star}$ are both
defined in (66). The idea is to let $\bm{\gamma}^{\star\tau}$ approximate the
singular values of $\bm{T}^{\star}\times_{3}\bm{\theta}^{\tau}$; this can be
seen, for instance, via the following calculation:

$\displaystyle\bm{T}^{\star}\times_{3}\bm{\theta}^{\tau}=\sum_{i=1}^{r}\big{\langle}\bm{u}_{i}^{\star},\bm{\theta}^{\tau}\big{\rangle}\bm{u}_{i}^{\star}\bm{u}_{i}^{\star\top}=\sum_{i=1}^{r}\underset{=\gamma_{i}^{\star\tau}}{\underbrace{\lambda_{i}^{\star}\big{\langle}\overline{\bm{u}}_{i}^{\star},\bm{\theta}^{\tau}\big{\rangle}}}\,\overline{\bm{u}}_{i}^{\star}\overline{\bm{u}}_{i}^{\star\top},$
(73)

where $\\{\overline{\bm{u}}_{i}^{\star}\\}_{i=1}^{r}$ — which are assumed to
be incoherent (or nearly orthogonal to each other) — can be approximately
viewed as the singular vectors of
$\bm{T}^{\star}\times_{3}\bm{\theta}^{\tau}$.

If we want our spectral estimate to be accurate, we would need to be assured
that the two largest entries of $\bm{\gamma}^{\star\tau}$ (in magnitude) are
sufficiently separated.

######  Lemma 5.12.

Instate the assumptions of Theorem 5.9. Define
$\Delta_{1}^{\tau}:=\gamma_{1}^{\star\tau}-\max_{1<i\leq
r}\big{|}\gamma_{i}^{\star\tau}\big{|}$ for each $1\leq\tau\leq L$ and let
$\Delta_{1}^{\left(1\right)}\geq\Delta_{1}^{\left(2\right)}\geq\dots\geq\Delta_{1}^{\left(L\right)}$
denote the order statistics of
$\big{\\{}\Delta_{1}^{\tau}\big{\\}}_{\tau=1}^{L}$ (in descending order). Fix
any arbitrary small constant $\delta>0$. With probability greater than
$1-\delta/r$, one has

$\displaystyle\Delta_{1}^{\left(1\right)}$
$\displaystyle\gtrsim\lambda_{\min}^{\star},$ (74a)
$\displaystyle\Delta_{1}^{\left(1\right)}-\Delta_{1}^{\left(2\right)}$
$\displaystyle\gtrsim\frac{\lambda_{\min}^{\star}}{r\sqrt{\log d}},$ (74b)

Additionally, for any fixed vector $\bm{v}\in\mathbb{R}^{r}$, with probability
at least $1-O\left(d^{-10}\right)$, for all $1\leq\tau\leq L$, one has

$\displaystyle\gamma_{1}^{\star\tau}$ $\displaystyle\lesssim\sqrt{\log
d}\,\lambda_{\max}^{\star},$ (75a)
$\displaystyle\big{\|}\bm{\gamma}^{\star\tau}\big{\|}_{2}$
$\displaystyle\lesssim\sqrt{r\log d}\,\lambda_{\max}^{\star},$ (75b)
$\displaystyle\big{|}\big{\langle}\bm{v},\bm{\gamma}^{\star\tau}\big{\rangle}\big{|}$
$\displaystyle\lesssim\left\|\bm{v}\right\|_{2}\sqrt{\log
d}\,\lambda_{\max}^{\star}.$ (75c)

###### Proof.

See Appendix B.1. ∎

Lemma 5.12 demonstrates that there exists some $\tau\in\left[L\right]$ such
that $\gamma_{1}^{\star\tau}-\max_{1<i\leq
r}\big{|}\gamma_{i}^{\star\tau}\big{|}\gtrsim\lambda_{\min}^{\star}$. This
means that $\overline{\bm{u}}_{1}^{\star}$ exhibits the largest correlation
with the random projection $\theta$, which further implies that
$\overline{\bm{u}}_{1}^{\star}$ is the largest singular vector of
$\bm{T}^{\star}\times_{3}\bm{\theta}^{\tau}$ with a considerable spectral gap
(as we will show shortly). With the desired spectral gap in place, we are
ready to look at the eigenvectors / singular vectors of interest. To this end,
we find it convenient to introduce another auxiliary vector
$\overline{\bm{u}}^{\tau}$, defined as the leading singular vector of
$\bm{M}^{\tau}$ (cf. (70c)) obeying

$\displaystyle\big{\langle}\overline{\bm{u}}^{\tau},\overline{\bm{u}}^{\star}_{1}\big{\rangle}\geq
0.$ (76)

The careful reader would immediately notice the similarity between
$\overline{\bm{u}}^{\tau}$ and $\bm{\nu}^{\tau}$ except for their global
signs; namely, we determine the global sign of $\overline{\bm{u}}^{\tau}$
based on the ground truth information (76), but pick the global sign for
$\bm{\nu}^{\tau}$ solely based on the observed data (cf. Algorithm 3).
Fortunately, the vectors $\overline{\bm{u}}^{\tau}$ and $\bm{\nu}^{\tau}$
provably coincide, namely,

$\overline{\bm{u}}^{\tau}=\bm{\nu}^{\tau},$ (77)

as we shall demonstrate momentarily in Lemma 5.16. In a similar way, we also
denote by $\overline{\bm{u}}^{\tau,\left(m\right)}$ the leading singular
vector of $\bm{M}^{\tau,\left(m\right)}$ defined in (70d) such that

$\displaystyle\big{\langle}\overline{\bm{u}}^{\tau,\left(m\right)},\overline{\bm{u}}^{\star}_{1}\big{\rangle}\geq
0.$ (78)

Lemma 5.16 also shows that
$\overline{\bm{u}}^{\tau,\left(m\right)}=\bm{\nu}^{\tau,\left(m\right)}$.

We shall now take a detour to look at $\overline{\bm{u}}^{\tau}$, which in
turn would help us understand $\bm{\nu}^{\tau}$. We shall first demonstrate
that $\overline{\bm{u}}^{\tau}$ (and hence $\bm{\nu}^{\tau}$) is sufficiently
close to the corresponding true factor in the $\ell_{2}$ sense.

######  Lemma 5.13.

Instate the assumptions of Theorem 5.9. Let $\overline{\bm{u}}^{\tau}$ and
$\overline{\bm{u}}_{1}^{\star}$ be as defined in (76) and (66), respectively.
Define $\mathcal{A}$ to be the event such that
$\gamma_{1}^{\star\tau}-\max_{1<i\leq
r}\big{|}\gamma_{i}^{\star\tau}\big{|}\gtrsim\lambda_{\min}^{\star}$ and the
condition (75) hold. Then conditional on this event $\mathcal{A}$, with
probability exceeding $1-O\left(d^{-11}\right)$ one has

$\displaystyle\big{\|}\overline{\bm{u}}^{\tau}-\overline{\bm{u}}_{1}^{\star}\big{\|}_{2}$
$\displaystyle\lesssim\underbrace{\frac{\mu r\log
d}{d\sqrt{p}}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{rd\log^{2}d}{p}}+\sqrt{\frac{\mu
r\log d}{d}}}_{=:\,\mathcal{E}_{\mathsf{proj}}}.$ (79)

###### Proof.

See Appendix B.2. ∎

Thus far, we have focused on the $\ell_{2}$ estimation errors. In order to
further quantify the $\ell_{\infty}$ estimation errors, we need to resort to
the leave-one-out estimates $\overline{\bm{u}}^{\tau,\left(m\right)}$ ($1\leq
m\leq d$). Specifically, we shall justify in the following two lemmas that:
(1) the $m$-th leave-one-out estimate
$\overline{\bm{u}}^{\tau,\left(m\right)}$ is close to the truth at least in
the $m$-th coordinate, and (2) the vector $\overline{\bm{u}}^{\tau}$ is
extremely close to each of the leave-one-out estimates
$\overline{\bm{u}}^{\tau,\left(m\right)}$ ($1\leq m\leq d$). These two
observations taken collectively translate to the desired entrywise error
control of $\overline{\bm{u}}^{\tau}$. Here, we recall that the global sign of
$\overline{\bm{u}}^{\tau,\left(m\right)}$ (cf. (78)) and the global sign of
$\overline{\bm{u}}^{\tau}$ (cf. (76)) are defined in a similar fashion, both
using the ground truth information.

######  Lemma 5.14.

Instate the assumptions of Theorem 5.9. Define $\mathcal{A}$ to be the event
such that $\gamma_{1}^{\star\tau}-\max_{1<i\leq
r}\big{|}\gamma_{i}^{\star\tau}\big{|}\gtrsim\lambda_{\min}^{\star}$ and the
condition (75) hold. Then conditional on this event $\mathcal{A}$, one has,
with probability exceeding $1-O\left(d^{-10}\right)$, that

$\displaystyle\Big{|}\big{[}\overline{\bm{u}}^{\tau,\left(m\right)}-\overline{\bm{u}}_{1}^{\star}\big{]}_{m}\Big{|}$
$\displaystyle\lesssim\mathcal{E}_{\mathsf{op}}\sqrt{\frac{\mu r\log d}{d}},$
(80)

holds for all $m\in\left[d\right]$, where $\mathcal{E}_{\mathsf{op}}$ is
defined as follows:

$\displaystyle\mathcal{E}_{\mathsf{op}}:=\frac{\sqrt{\mu
r}\,\log^{3}d}{d^{3/2}p}+\frac{\mu
r\log^{5/2}d}{d\sqrt{p}}+\frac{\sigma}{\lambda_{\min}^{\star}}\frac{\log^{7/2}d}{p}+\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{rd\log^{5}d}{p}}+\sqrt{\frac{\mu
r\log d}{d}}.$ (81)

###### Proof.

See Appendix B.4. ∎

######  Lemma 5.15.

Instate the assumptions of Theorem 5.9. Define $\mathcal{A}$ to be the event
such that $\gamma_{1}^{\star\tau}-\max_{1<i\leq
r}\big{|}\gamma_{i}^{\star\tau}\big{|}\gtrsim\lambda_{\min}^{\star}$ and the
condition (75) hold. Then conditional on this event $\mathcal{A}$, one has,
with probability at least $1-O\left(d^{-10}\right)$, for all
$m\in\left[d\right]$:

$\displaystyle\big{\|}\overline{\bm{u}}^{\tau}-\overline{\bm{u}}^{\tau,\left(m\right)}\big{\|}_{2}$
$\displaystyle\lesssim\mathcal{E}_{\mathsf{loo}}\sqrt{\frac{\mu r\log d}{d}},$
(82)
$\displaystyle\big{\|}\overline{\bm{u}}^{\tau}-\overline{\bm{u}}_{1}^{\star}\big{\|}_{\infty}$
$\displaystyle\lesssim\left(\mathcal{E}_{\mathsf{op}}+\mathcal{E}_{\mathsf{loo}}\right)\sqrt{\frac{\mu
r\log d}{d}},$ (83)
$\displaystyle\big{|}\lambda_{\tau}-\lambda_{\tau}^{\left(m\right)}\big{|}$
$\displaystyle\lesssim\mathcal{E}_{\mathsf{loo}}\sqrt{\frac{\mu r\log
d}{d}}\,\lambda_{\max}^{\star},$ (84)

where $\mathcal{E}_{\mathsf{op}}$ and $\mathcal{E}_{\mathsf{loo}}$ are defined
in (81) and (65), respectively.

###### Proof.

See Appendix B.6. ∎

Next, we turn to the estimation accuracy regarding the size of the tensor
factors and show that $\lambda_{\tau}$ (produced in Algorithm 3) is close to
the truth as well. As it turns out, a byproduct of this step reveals that
$\bm{\nu}^{\tau}=\overline{\bm{u}}^{\tau}$ and
$\bm{\nu}^{\tau,\left(m\right)}=\overline{\bm{u}}^{\tau,\left(m\right)}$,
where $\overline{\bm{u}}^{\tau}$ and $\overline{\bm{u}}^{\tau,\left(m\right)}$
are an auxiliary vectors defined in (76) and (78), respectively.

######  Lemma 5.16.

Instate the assumptions of Theorem 5.9. Assume that the results in Lemma 5.13,
Lemma 5.14 and Lemma 5.15 hold. Then with probability at least
$1-O\left(d^{-10}\right)$, one has

$\displaystyle\big{|}\lambda_{\tau}-\lambda_{1}^{\star}\big{|}\lesssim\mathcal{E}_{\mathsf{proj}}\lambda_{1}^{\star}.$
(85)

In particular, one has

$\displaystyle\bm{\nu}^{\tau}=\overline{\bm{u}}^{\tau}\qquad\text{and}\qquad\bm{\nu}^{\tau,\left(m\right)}=\overline{\bm{u}}^{\tau,\left(m\right)},\quad
1\leq m\leq d.$ (86)

###### Proof.

See Appendix B.8. ∎

Thus far, we have only proved that one can find a reliable estimate for each
tensor factor within $L$ random trials, provided that $L$ is sufficiently
large. To finish up, it remains to show that the pruning procedure
$\textsc{Prune}($) is capable of returning a rough estimate for each tensor
factor without duplication. This is accomplished in the following lemma.

######  Lemma 5.17.

Instate the assumptions of Theorem 5.9. On the event that the results in Lemma
5.13, Lemma 5.14, Lemma 5.15 and Lemma 5.16 hold for all $1\leq i\leq r$,
there exists a permutation $\pi(\cdot):[d]\mapsto[d]$ such that: for each
$1\leq i\leq r$, $\big{(}\lambda_{i},\bm{w}^{i}\big{)}$ and
$\big{(}\lambda_{\pi(i)}^{\star},\overline{\bm{u}}^{\star}_{\pi(i)}\big{)}$
satisfy (68a), (68b) and (68c); $\big{(}\lambda_{i},\bm{w}^{i}\big{)}$ and
$\big{(}\lambda_{i}^{\left(m\right)},\bm{w}^{i,\left(m\right)}\big{)}_{i=1}^{r}$
obey (69a), (69b) and (69c) for all $1\leq m\leq d$, where
$\big{\\{}\lambda_{i},\bm{w}^{i}\big{\\}}_{i=1}^{r}$ and
$\big{\\{}\lambda_{i}^{\left(m\right)},\bm{w}^{i,\left(m\right)}\big{\\}}_{i=1}^{r}$
are outputs of Algorithm 3 and Algorithm 6, respectively.

###### Proof.

See Appendix B.9. ∎

##  6 Discussion

The current paper uncovers the possibility of efficiently and stably
completing a low-CP-rank tensor from partial and noisy entries. Perhaps
somewhat unexpectedly, despite the high degree of nonconvexity, this problem
can be solved to optimal statistical accuracy within nearly linear time,
provided that the tensor of interest is well-conditioned, incoherent, and of
constant rank. To the best of our knowledge, this intriguing message has not
been shown in the prior literature.

Moving forward, one pressing issue is to understand how to improve the
algorithmic and theoretical dependency upon the tensor rank $r$ of the
proposed method. Ideally one would desire a fast algorithm whose sample
complexity scales as $rd^{1.5}$, an order that is provably achievable by the
sum-of-squares hierarchy. Additionally, in contrast to the matrix counterpart
where the rank is upper bounded by the matrix dimension, the tensor CP rank is
allowed to rise above $d$, which is commonly referred to as the over-complete
case. Unfortunately, our current initialization scheme (i.e. the spectral
method) fails to work unless $r<d$, and our local analysis for GD falls of
accommodating the scenario with $r>d$. It would be of great interest to
develop more powerful algorithms — in addition to more refined analysis — to
tackle such an important over-complete regime.

Another tantalizing research direction is the exploration of landscape design
for tensor completion. As our heuristic discussions as well as other prior
work (e.g. [RM14]) suggest, randomly initialized gradient descent tailored to
(4) seems unlikely to work, unless the sample size is significantly larger
than the computational limit. This might mean either that there exist spurious
local minima in the natural nonconvex least squares formulation (4), or that
the optimization landscape of (4) is too flat around some saddle points and
hence not amenable to fast computation. It would be interesting to investigate
what families of loss functions allow us to rule out bad local minima and
eliminate the need of careful initialization, which might be better suited for
tensor recovery problems.

Finally, in statistical inference and decision making, one might not be simply
satisfied with obtaining a reliable estimate for each missing entry, but would
also like to report a short confidence interval which is likely to contain the
true entry. This boils down to the fundamental task of uncertainty
quantification for tensor completion, which we leave to future investigation.

## Acknowledgements

Y. Chen is supported in part by the AFOSR YIP award FA9550-19-1-0030, by the
ONR grant N00014-19-1-2120, by the ARO grants W911NF-20-1-0097 and
W911NF-18-1-0303, by the NSF grants CCF-1907661, IIS-1900140 and DMS-2014279,
and by the Princeton SEAS innovation award. H. V. Poor is supported in part by
the NSF grant DMS-1736417. C. Cai is supported in part by Gordon Y. S. Wu
Fellowships in Engineering. This work was done in part while Y. Chen was
visiting the Kavli Institute for Theoretical Physics (supported in part by NSF
grant PHY-1748958). We thank Lanqing Yu for many helpful discussions, and
thank Yuling Yan for proofreading the paper.

##  Appendix A Proofs for local convergence of GD

In this section, we establish the key lemmas concerning the convergence
properties of GD. As one can easily see, treating $\\{E_{i,j,k}\\}_{1\leq
i,j,k\leq d}$ (resp. $\\{\chi_{i,j,k}\\}_{1\leq i,j,k\leq d}$) as independent
random variables — which leads to asymmetric versions of $\bm{E}$ and $\Omega$
— does not affect the order of our results at all. In light of this, we shall
adopt such an independent assumption whenever it simplifies our presentation.

###  A.1 Proof of Lemma 5.1

For notational convenience, for any matrix
$\bm{M}=\left[\bm{m}_{1},\dots,\bm{m}_{r}\right]\in\mathbb{R}^{d\times r}$,
let

$\displaystyle\widetilde{\bm{M}}:=\left[\bm{m}_{1}\otimes\bm{m}_{1},\dots,\bm{m}_{r}\otimes\bm{m}_{r}\right]\in\mathbb{R}^{d^{2}\times
r},$ (87)

where for any $\bm{a},\bm{b}\in\mathbb{R}^{d}$ we denote
$\bm{a}\otimes\bm{b}:={\footnotesize\left[\begin{array}[]{c}a_{1}\bm{b}\\\
\vdots\\\ a_{d}\bm{b}\end{array}\right]}\in\mathbb{R}^{d^{2}}$.

From the Hessian expression (41), one can decompose

$\displaystyle\mathsf{vec}(\bm{V})^{\top}\nabla^{2}f(\bm{U})\,\mathsf{vec}(\bm{V})$
$\displaystyle=\underbrace{\frac{1}{3p}\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\big{(}\bm{v}_{s}\otimes\bm{u}_{s}^{\otimes
2}+\bm{u}_{s}\otimes\bm{v}_{s}\otimes\bm{u}_{s}+\bm{u}_{s}^{\otimes
2}\otimes\bm{v}_{s}\big{)}\Big{)}\Big{\|}_{\mathrm{F}}^{2}-\frac{1}{3p}\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\big{(}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}+\bm{u}_{s}^{\star}\otimes\bm{v}_{s}\otimes\bm{u}_{s}^{\star}+\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\big{)}\Big{)}\Big{\|}_{\mathrm{F}}^{2}}_{=:\,\alpha_{1}}$
$\displaystyle\quad+\underbrace{\frac{1}{3p}\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\big{(}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}+\bm{u}_{s}^{\star}\otimes\bm{v}_{s}\otimes\bm{u}_{s}^{\star}+\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\big{)}\Big{)}\Big{\|}_{\mathrm{F}}^{2}-\frac{1}{3}\,\Big{\|}\sum_{s\in[r]}\big{(}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}+\bm{u}_{s}^{\star}\otimes\bm{v}_{s}\otimes\bm{u}_{s}^{\star}+\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\big{)}\Big{\|}_{\mathrm{F}}^{2}}_{=:\,\alpha_{2}}$
$\displaystyle\quad+2\underbrace{\Big{\langle}\frac{1}{p}\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{s\in[r]}\bm{u}_{s}^{\otimes
3}-\bm{T}^{\star}\Big{)},\sum\nolimits_{s\in[r]}\bm{v}_{s}^{\otimes
2}\otimes\bm{u}_{s}\Big{\rangle}}_{=:\,\alpha_{3}}+\underbrace{\frac{1}{3}\,\Big{\|}\sum_{s\in[r]}\big{(}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}+\bm{u}_{s}^{\star}\otimes\bm{v}_{s}\otimes\bm{u}_{s}^{\star}+\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\big{)}\Big{\|}_{\mathrm{F}}^{2}}_{=:\,\alpha_{4}}.$

In what follows, we shall bound each of the above terms separately.

####  A.1.1 Bounding $\alpha_{4}$

With regards to $\alpha_{4}$, by symmetry we have

$\displaystyle\alpha_{4}=\Big{\|}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{\|}_{\mathrm{F}}^{2}+2\,\Big{\langle}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2},\sum\nolimits_{s\in[r]}\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\Big{\rangle}.$ (88)

In order to control (88), we first see that

$\Big{\|}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{\|}_{\mathrm{F}}=\big{\|}\bm{V}\widetilde{\bm{U}}^{\star\top}\big{\|}_{\mathrm{F}},$
(89)

where $\widetilde{\bm{U}}^{\star}$ is as defined in (87). Similar to the proof
of Lemma D.1, we can use the fact that $\langle\bm{u}_{i}^{\star\otimes
2},\bm{u}_{j}^{\star\otimes
2}\rangle=\langle\bm{u}_{i}^{\star},\bm{u}_{j}^{\star}\rangle^{2}$ and (8c) to
deduce that

$\displaystyle\sigma_{\min}\big{(}\widetilde{\bm{U}}^{\star}\big{)}=\lambda_{\min}^{\star
2/3}\left(1+o\left(1\right)\right)\qquad\text{and}\qquad\mathrm{\sigma}_{\max}\big{(}\widetilde{\bm{U}}^{\star}\big{)}=\lambda_{\max}^{\star
2/3}\left(1+o\left(1\right)\right),$ (90)

provided that $r\ll d/\mu$. This implies that

$\displaystyle\frac{19}{20}\lambda_{\min}^{\star
2/3}\left\|\bm{V}\right\|_{\mathrm{F}}\leq\sigma_{\min}\big{(}\widetilde{\bm{U}}^{\star}\big{)}\left\|\bm{V}\right\|_{\mathrm{F}}\leq\big{\|}\bm{V}\widetilde{\bm{U}}^{\star\top}\big{\|}_{\mathrm{F}}\leq\sigma_{\max}\big{(}\widetilde{\bm{U}}^{\star}\big{)}\left\|\bm{V}\right\|_{\mathrm{F}}\leq\frac{11}{10}\lambda_{\max}^{\star
2/3}\left\|\bm{V}\right\|_{\mathrm{F}}.$ (91)

(1) Speaking of an upper bound on $\alpha_{4}$, we can invoke the Cauchy-
Schwarz inequality followed by (91) to reach

$\displaystyle\alpha_{4}\leq
3\,\Big{\|}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{\|}_{\mathrm{F}}^{2}=3\,\big{\|}\bm{V}\widetilde{\bm{U}}^{\star\top}\big{\|}_{\mathrm{F}}^{2}\leq\frac{7}{2}\lambda_{\max}^{\star
4/3}\left\|\bm{V}\right\|_{\mathrm{F}}^{2}.$ (92)

(2) When it comes to lower bounding $\alpha_{4}$, the main step boils down to
controlling the inner product term in (88). Applying the Cauchy-Schwartz
inequality gives that

$\displaystyle\Big{\langle}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2},\sum\nolimits_{s\in[r]}\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\Big{\rangle}$
$\displaystyle=\sum_{s\in[r]}\left\langle\bm{v}_{s},\bm{u}_{s}^{\star}\right\rangle^{2}\left\|\bm{u}_{s}^{\star}\right\|_{2}^{2}+\sum_{s_{1}\neq
s_{2}}\left\langle\bm{v}_{s_{1}},\bm{u}^{\star}_{s_{2}}\right\rangle\left\langle\bm{u}^{\star}_{s_{1}},\bm{v}_{s_{2}}\right\rangle\left\langle\bm{u}^{\star}_{s_{1}},\bm{u}^{\star}_{s_{2}}\right\rangle$
$\displaystyle\geq-\max_{s_{1}\neq
s_{2}}\left|\left\langle\bm{u}^{\star}_{s_{1}},\bm{u}^{\star}_{s_{2}}\right\rangle\right|\sum_{s_{1}\neq
s_{2}}\left\|\bm{v}_{s_{1}}\right\|_{2}\left\|\bm{u}^{\star}_{s_{1}}\right\|_{2}\left\|\bm{v}_{s_{2}}\right\|_{2}\left\|\bm{u}^{\star}_{s_{2}}\right\|_{2}$
$\displaystyle\geq-\max_{s_{1}\neq
s_{2}}\left|\left\langle\bm{u}^{\star}_{s_{1}},\bm{u}^{\star}_{s_{2}}\right\rangle\right|\Big{(}\sum\nolimits_{s\in[r]}\left\|\bm{v}_{s}\right\|_{2}\left\|\bm{u}^{\star}_{s}\right\|_{2}\Big{)}^{2}$
$\displaystyle\overset{\text{(i)}}{\geq}-\max_{s_{1}\neq
s_{2}}\left|\left\langle\bm{u}^{\star}_{s_{1}},\bm{u}^{\star}_{s_{2}}\right\rangle\right|\left\|\bm{U}^{\star}\right\|_{\mathrm{F}}^{2}\left\|\bm{V}\right\|_{\mathrm{F}}^{2}$
$\displaystyle\geq-r\sqrt{\frac{\mu}{d}}\,\lambda_{\max}^{\star
4/3}\left\|\bm{V}\right\|_{\mathrm{F}}^{2}\geq-\frac{1}{40}\lambda_{\min}^{\star
4/3}\left\|\bm{V}\right\|_{\mathrm{F}}^{2},$

where (i) comes from Cauchy-Schwarz, and the last line follows from (8c), (8d)
as well as the condition that $r\ll\sqrt{d/\mu}$ and $\kappa\asymp 1$.
Therefore, we can lower bound $\alpha_{4}$ by (with the assistance of (91))

$\displaystyle\alpha_{4}$
$\displaystyle=\,\Big{\|}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{\|}_{\mathrm{F}}^{2}+2\,\Big{\langle}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2},\sum_{s\in[r]}\bm{u}_{s}^{\star\otimes 2}\otimes\bm{v}_{s}\Big{\rangle}$
$\displaystyle\geq\frac{19}{20}\lambda_{\min}^{\star
4/3}\left\|\bm{V}\right\|_{\mathrm{F}}^{2}-\frac{1}{20}\lambda_{\min}^{\star
4/3}\left\|\bm{V}\right\|_{\mathrm{F}}^{2}\leavevmode\nobreak\
\geq\leavevmode\nobreak\ \frac{9}{10}\lambda_{\min}^{\star
4/3}\left\|\bm{V}\right\|_{\mathrm{F}}^{2}.$ (93)

####  A.1.2 Bounding $\alpha_{1}$

When it comes to $\alpha_{1}$, we can expand

$\displaystyle\Big{\|}$
$\displaystyle\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\otimes
2}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{s\in[r]}\bm{u}_{s}\otimes\bm{v}_{s}\otimes\bm{u}_{s}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{s\in[r]}\bm{u}_{s}^{\otimes
2}\otimes\bm{v}_{s}\Big{)}\Big{\|}_{\mathrm{F}}^{2}$
$\displaystyle=3\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\otimes
2}\Big{)}\Big{\|}_{\mathrm{F}}^{2}+6\,\Big{\langle}\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\otimes
2}\Big{)},\mathcal{P}_{\Omega}\Big{(}\sum\nolimits_{s\in[r]}\bm{u}_{s}^{\otimes
2}\otimes\bm{v}_{s}\Big{)}\Big{\rangle};$

we can decompose
$\big{\|}\mathcal{P}_{\Omega}\big{(}\sum\nolimits_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\big{)}+\mathcal{P}_{\Omega}\big{(}\sum\nolimits_{s\in[r]}\bm{u}_{s}^{\star}\otimes\bm{v}_{s}\otimes\bm{u}_{s}^{\star}\big{)}+\mathcal{P}_{\Omega}\big{(}\sum\nolimits_{s\in[r]}\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\big{)}\big{\|}_{\mathrm{F}}^{2}$ in a similar way. As a
consequence,

$\displaystyle\alpha_{1}$
$\displaystyle=\underbrace{\frac{2}{p}\Bigg{(}\Big{\langle}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\otimes
2}\Big{)},\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\otimes
2}\otimes\bm{v}_{s}\Big{)}\Big{\rangle}-\Big{\langle}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{)},\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\Big{)}\Big{\rangle}\Bigg{)}}_{=:\,\beta_{1}}$
$\displaystyle\quad+\underbrace{\frac{1}{p}\Bigg{(}\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\otimes
2}\Big{)}\Big{\|}_{\mathrm{F}}^{2}-\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{)}\Big{\|}_{\mathrm{F}}^{2}\Bigg{)}}_{=:\,\beta_{2}}.$

We will derive an upper bound on $\beta_{1}$ in the sequel; the same method
immediately applies to $\beta_{2}$.

For notational convenience, let us define

$\displaystyle\bm{\Delta}:=\bm{U}-\bm{U}^{\star},\quad\bm{\Delta}_{s}:=\bm{u}_{s}-\bm{u}_{s}^{\star},\quad\widetilde{\bm{\Delta}}:=[\bm{\Delta}_{1}\otimes\bm{\Delta}_{1},\cdots,\bm{\Delta}_{r}\otimes\bm{\Delta}_{r}]\in\mathbb{R}^{d^{2}\times
r}.$ (94)

Then one can write

$\displaystyle\frac{1}{2}\beta_{1}$
$\displaystyle=\frac{1}{p}\,\Big{\langle}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes(\bm{\Delta}_{s}+\bm{u}_{s}^{\star})^{\otimes
2}\Big{)},\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}(\bm{\Delta}_{s}+\bm{u}_{s}^{\star})^{\otimes
2}\otimes\bm{v}_{s}\Big{)}\Big{\rangle}-\frac{1}{p}\,\Big{\langle}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{)},\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\Big{)}\Big{\rangle}$
$\displaystyle=\frac{1}{p}\,\Big{\langle}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{)},\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\star}\otimes\bm{\Delta}_{s}\otimes\bm{v}_{s}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{\Delta}_{s}\otimes\bm{u}_{s}^{\star}\otimes\bm{v}_{s}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{\Delta}_{s}\otimes\bm{\Delta}_{s}\otimes\bm{v}_{s}\Big{)}\Big{\rangle}$
$\displaystyle\quad+\frac{1}{p}\,\Big{\langle}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star}\otimes\bm{\Delta}_{s}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{\Delta}_{s}\otimes\bm{u}_{s}^{\star}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{\Delta}^{\otimes
2}_{s}\Big{)},$
$\displaystyle\qquad\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\star}\otimes\bm{\Delta}_{s}\otimes\bm{v}_{s}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{\Delta}_{s}\otimes\bm{u}_{s}^{\star}\otimes\bm{v}_{s}\Big{)}+\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{\Delta}_{s}\otimes\bm{\Delta}_{s}\otimes\bm{v}_{s}\Big{)}\Big{\rangle}.$

Apply the Cauchy-Schwartz inequality to yield that

$\displaystyle\left|\beta_{1}\right|$
$\displaystyle\lesssim\frac{1}{p}\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star\otimes
2}\Big{)}\Big{\|}_{\mathrm{F}}\Big{(}2\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\star}\otimes\bm{\Delta}_{s}\otimes\bm{v}_{s}\Big{)}\Big{\|}_{\mathrm{F}}+\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{\Delta}_{s}\otimes\bm{\Delta}_{s}\otimes\bm{v}_{s}\Big{)}\Big{\|}_{\mathrm{F}}\Big{)}$
$\displaystyle\quad+\frac{1}{p}\Big{(}\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\star\otimes
2}\otimes\bm{v}_{s}\Big{)}\Big{\|}_{\mathrm{F}}+2\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{u}_{s}^{\star}\otimes\bm{\Delta}_{s}\otimes\bm{v}_{s}\Big{)}\Big{\|}_{\mathrm{F}}+\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{\Delta}_{s}\otimes\bm{\Delta}_{s}\otimes\bm{v}_{s}\Big{)}\Big{\|}_{\mathrm{F}}\Big{)}$
$\displaystyle\quad\quad\cdot\Big{(}2\,\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{u}_{s}^{\star}\otimes\bm{\Delta}_{s}\Big{)}\Big{\|}_{\mathrm{F}}+\Big{\|}\mathcal{P}_{\Omega}\Big{(}\sum_{s\in[r]}\bm{v}_{s}\otimes\bm{\Delta}^{\otimes
2}_{s}\Big{)}\Big{\|}_{\mathrm{F}}\Big{)}.$ (95)

Before we bound the above quantities, we pause to make the following
observations. In view of the assumptions of this lemma that $\delta\ll
1/\sqrt{r}\leq 1$, the following holds for all $i\in\left[r\right]$:

$\displaystyle\left\|\bm{\Delta}_{i}\right\|_{2}$
$\displaystyle\leq\left\|\bm{U}-\bm{U}^{\star}\right\|_{\mathrm{F}}\leq\delta\left\|\bm{U}^{\star}\right\|_{\mathrm{F}}\leq\delta\sqrt{r}\,\lambda_{\max}^{\star
1/3}\ll\lambda_{\max}^{\star 1/3},$ (96a)
$\displaystyle\left\|\bm{\Delta}_{i}\right\|_{\infty}$
$\displaystyle\leq\left\|\bm{U}-\bm{U}^{\star}\right\|_{2,\infty}\leq\delta\left\|\bm{U}^{\star}\right\|_{2,\infty}\leq\delta\sqrt{\frac{\mu
r}{d}}\,\lambda_{\max}^{\star 1/3}\ll\sqrt{\frac{\mu}{d}}\lambda_{\max}^{\star
1/3},$ (96b) $\displaystyle\left\|\bm{u}_{i}\right\|_{2}$
$\displaystyle\leq\left\|\bm{u}_{i}^{\star}\right\|_{2}+\left\|\bm{\Delta}_{i}\right\|_{2}\leq
2\lambda_{\max}^{\star 1/3},$ (96c)
$\displaystyle\left\|\bm{u}_{i}\right\|_{\infty}$
$\displaystyle\leq\left\|\bm{u}_{i}^{\star}\right\|_{\infty}+\left\|\bm{\Delta}_{i}\right\|_{\infty}\leq
2\sqrt{\frac{\mu}{d}}\,\lambda_{\max}^{\star 1/3},$ (96d)

Consequently, we also know that

$\displaystyle\big{\|}\widetilde{\bm{\Delta}}\big{\|}_{2,\infty}\leq\max_{1\leq
i\leq
r}\left\|\bm{\Delta}_{i}\right\|_{\infty}\left\|\bm{\Delta}\right\|_{2,\infty}\leq\delta^{2}\left\|\bm{U}^{\star}\right\|_{2,\infty}^{2}.$
(97)

Now, we proceed to prove the claim. Let us define
$\mathcal{S}_{i}:=\\{\underline{j}\in[d]^{2}\mid\chi_{ij_{1}j_{2}}=1\\}$ for
each $i\in[d]$. Applying the Chernoff bound and the union bound yields that:
with probability at least $1-O(d^{-10})$ one has

$\displaystyle\max_{i\in\left[d\right]}|\mathcal{S}_{i}|\lesssim d^{2}p.$ (98)

provided $p\gg d^{-2}\log d$. It then follows from the Cauchy-Schwarz
inequality that