Nonconvex Low-Rank Tensor Completion from Noisy Data

Introduction and motivation

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.

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.

Algorithm and main results

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

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

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.

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

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.

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:

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

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.

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.

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

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]).

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:

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

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$.

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

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$,

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.

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

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.

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

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.

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.

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.

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)

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.

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.

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.

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.

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

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$,

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.

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

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.

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

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.

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}_{},\cdots,\bm{U}^{0}_{[t_{\mathsf{init}}]}$); select the one achieving the smallest empirical loss, namely

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

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$.

4 Notation

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

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$.

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.

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.

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

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

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,

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

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.

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:

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

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 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

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

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

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

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))

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

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.

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,

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.

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.

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.

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

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

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

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.

Suppose that the sample complexity and the rank satisfy

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

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

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.

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:

$\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:

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.

1.4 Key lemmas

Key hypotheses for the gradient update stage:

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.

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

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.

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})$,

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.

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

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.

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

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.

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

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.

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

2 Analysis for initialization: Part 1 (subspace estimation)

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}}$.

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

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

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

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.

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$:

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

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

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}$.

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

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

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

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.

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

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

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

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

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.

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$:

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.

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$.

3.2 Analysis

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

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

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

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.

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}$.

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.

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:

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.

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

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

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,

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

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

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

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

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

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]$:

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

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.

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

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.

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.

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.

From the Hessian expression (41), one can decompose

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

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

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

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

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

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

(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

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))

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

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

Apply the Cauchy-Schwartz inequality to yield that

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]$:

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

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

where (i) arises from (97) and (98). In a similar manner, we can derive

under the sample size assumptin that $p\gg\mu^{2}r^{2}d^{-2}\log d$. Here the last inequality makes use of (90). It is self-evident that the above bounds also hold for quantities that appear in (95). Since $0<\delta\ll 1/\sqrt{r}<1$, we obtain

The same upper bound holds for any other $\beta_{2}$. Therefore, as long as $0<\delta\ll 1/(\mu\sqrt{r})$ and $\kappa\asymp 1$, we have

Regarding $\alpha_{2}$, applying [YZ16, Lemma 5] (with slight modification, which we omit here for brevity) implies that: if $p\gg\mu^{2}r^{2}d^{-2}\log d$, then with probability exceeding $1-O(d^{-10})$,

Here the last inequality arises from (92).

We now move on to bounding $\alpha_{3}$. The triangle inequality gives

Recall the definitions of $\bm{\Delta}$ and $\bm{\Delta}_{i}$ in (94). Fix an arbitrary $s\in\left[r\right]$. From the definition of the operator norm and the triangle inequality, we can derive

With this lemma in mind, we are ready to derive that

Given that $p\gtrsim d^{-3/2}\log^{3}d$, applying Lemma D.2 indicates that with probability at least $1-O(d^{-10})$,

Moreover, it is straightforward to see that $\left\|\bm{1}^{\otimes 3}\right\|=\left\|\bm{1}\right\|_{2}^{3}=d^{3/2}$. Therefore, one has

Next, we turn to $\left\|\cdot\right\|_{\infty}$. We first expand

By symmetry, it suffices to control \big{\|}\sum_{i\in[r]}\bm{u}^{\star\otimes 2}_{i}\otimes\bm{\Delta}_{i}\big{\|}, \big{\|}\sum_{i\in[r]}\bm{\Delta}_{i}^{\otimes 2}\otimes\bm{u}^{\star}_{i}\big{\|} and \big{\|}\sum_{i\in[r]}\bm{\Delta}_{i}^{\otimes 3}\big{\|}. Let us look at the first term. Towards this, for each $\left(i,j,k\right)\in\left[d\right]^{3}$, we can use the Cauchy-Schwartz inequality to control

In a similar manner, we can control the remaining two terms by

Recall that $0<\delta\ll 1/r\leq 1$. Putting these together reveals that

where we use (96c) in the last step. In view of the condition that $\delta\ll 1/(\mu^{3/2}r)$ and the assumption $\kappa\asymp 1$, one has with probability greater than $1-O(d^{-10})$,

A.1.5 Putting all this together

Note that the above bounds hold uniformly for all $\bm{V}$. Therefore, combining upper bounds for $\alpha_{i}$ and the union bound, we conclude that with probability exceeding $1-O(d^{-10})$,

A.2 Proof of Lemma 5.2

From (47a) and (47c), we use the triangle inequality to obtain

Similarly, we can combine (47b) and (47c) to obtain that

A.3 Proof of Lemma 5.3

which motivates us to bound $\alpha_{1}$ and $\alpha_{2}$ separately.

(1) We start with $\alpha_{1}$, towards which we find it helpful to define

Given that $\nabla f_{\mathsf{clean}}\left(\bm{U}^{\star}\right)=\bm{0}$ (since $\bm{U}^{\star}$ is a global optimizer of $f_{\mathsf{clean}}$), we can use the fundamental theorem of calculus to obtain

From the hypothesis (47b) as well as our conditions that $\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log d}{p}}+\mathcal{E}_{\mathsf{local}}\ll\frac{1}{\mu^{3/2}r}$, we know that $\bm{U}^{t}(\tau)$ ($0\leq\tau\leq 1$) satisfies the conditions required in Lemma 5.1. Therefore, applying Lemma 5.1 gives that

Substitution into (116) indicates that: if 0<\eta\leq\lambda_{\min}^{\star 4/3}/\big{(}32\lambda_{\max}^{\star 8/3}\big{)}, then

which implies that (since $1-a/2\geq\sqrt{1-a}$ for $0<a<1$)

(2) We now turn to $\alpha_{2}$. To simplify presentation, we shall assume that $\{E_{i,j,k}\}_{i,j,k\in[d]}$ (resp. $\{\chi_{i,j,k}\}_{i,j,k\in[d]}$) are independent random variables. Fix an arbitrary $s\in[r]$ and $m\in[d]$. The $m$-th entry of $\mathcal{P}_{\Omega}\left(\bm{E}\right)\times_{1}\bm{u}^{t}_{s}\times_{2}\bm{u}^{t}_{s}$ can be expanded as follows:

Before continuing, we make the following observations: from the hypotheses (47a), (47b), (47c), as well as our assumption that $\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log d}{p}}+\mathcal{E}_{\mathsf{local}}\ll\frac{1}{\mu^{3/2}r}$, we see that the following holds for all $s\in\left[r\right]$:

With these estimates in place, we can upper bound the above four terms in (118) separately.

For $\beta_{1}$, we note that, by construction, $\bm{u}^{t,\left(m\right)}$ is independent of the $m$-th mode-3 slice of $\mathcal{P}_{\Omega}\left(\bm{E}\right)$. This tells us that

can be viewed as a sum of independent zero-mean random variables (conditional on $\mathcal{P}_{\Omega_{-m}}\left(\bm{E}\right)$). It is straightforward to compute that

where $\|\cdot\|_{\psi_{1}}$ denotes the sub-exponential norm and we use (119) in the last inequality. Applying the matrix Bernstein inequality [Kol11, Corollary 2.1] yields that: with probability $1-O(d^{-20})$,

For $\beta_{2}$, we first invoke [CW15, lemma 11] to demonstrate that: with probability at least $1-O(d^{-11})$,

where the last inequality follows from the noise condition. Clearly, the above bound holds for $\beta_{3}$ as well.

Taking together the above bounds and substituting them into (118), we obtain

Recognizing that this holds for any $m\in[d]$ and $s\in[r]$, one can sum over $m$ and $s$ to deduce that

(3) Combining (117) and (122) yields that: with probability at least $1-O(d^{-10})$,

with the proviso 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 that $C_{2}$ is sufficiently large.

A.4 Proof of Lemma 5.4

Fix an arbitrary $m\in[r]$. From the definition of $f^{\left(m\right)}$ in (45) and (46), we can show that

Before proceeding to bounding these terms, we pause to define

From (119) and hypothesis (47d), one can applying a similar argument as in (107) to find that

Now we begin to bound the terms in (123) separately.

(1) We start with $\alpha_{1}$. For any $0\leq\tau\leq 1$, define

The fundamental theorem of calculus yields

provided that 0<\eta\leq\lambda_{\min}^{\star 4/3}/\big{(}32\lambda_{\max}^{\star 8/3}\big{)}.

In what follows, we shall assume that $\{E_{i,j,k}\}_{i,j,k\in[d]}$ (resp. $\{\chi_{i,j,k}\}_{i,j,k\in[d]}$) are independent random variables in order to simplify presentation.

(2) The next step is to bound $\alpha_{2}$. For notational simplicity, define

In order to control the Frobenius norm of $\bm{V}^{t,\left(m\right)}$, we shall start by considering the $m$-th row of $\bm{V}^{t,\left(m\right)}$. In view of the definitions of $\mathcal{P}_{\Omega_{m}}$ and $\mathcal{P}_{m}$ (cf. Appendix 5.1.3), we can expand

We recognize that \big{\{}\bm{u}_{s}^{t,\left(m\right)}\big{\}}_{s=1}^{r} is independent of $\Omega_{m}$, making it convenient for us to upper bound \big{\|}\bm{V}^{t,\left(m\right)}_{m,:}\big{\|}_{2}. Specifically, for any $i,j\in[d]$, from (119) and (125) we have

where the last inequality holds due to (119) and (125). We then apply the matrix Bernstein inequality to yield that: with probability exceeding $1-O(d^{-20})$,

where the last step holds as long as $p\gg\mu^{2}d^{-2}\log d$.

Next, we turn to the $k$-th row of $\bm{V}^{t,\left(m\right)}$ for any $k\neq m$. For each $s\in[r]$, we have

Similar to the proof of Lemma D.9, we can show that with probability at least $1-O(d^{-20})$,

where the last inequality follows from (119), and (125). It is easily seen that the bound also holds for the summation over $i\neq m$. We can then use Cauchy-Schwartz to arrive at

Combining (128) with (127) and invoking the union bound, we conclude that: with probability exceeding $1-O(d^{-20})$,

for some absolute constant $C>0$, where we use the assumption that $\kappa\asymp 1$.

(3) For $\alpha_{3}$, following a similar argument for $\alpha_{2}$, we define

for each $s\in[r]$. The $m$-th row of $\bm{W}^{t,\left(m\right)}$ is a sum of independent zero-mean random vectors:

With (119) in place, it is easy to verify that

where the last inequality holds as long as $p\gg\mu d^{-2}\log^{3}d$.

for each $s\in[r]$. Arguing similarly as in the proof of Lemma D.10, we have with probability at least $1-O(d^{-20})$,

where the last inequality follows from (119). Additionally, the summation over $\{i:i\neq m\}$ can be controlled using the same argument. Therefore, we use Cauchy-Schwarz to find that

Combined with (130) and the assumption that $\kappa\asymp 1$, we obtain that

for some absolute constant $\widetilde{C}>0$.

(4) Regarding $\alpha_{4}$, we use the triangle inequality to show that for each $s\in[r]$,

where the last line follows from (119). From Corollary D.3, we can further upper bound

As a result, we sum over $s\in[r]$ and use the Cauchy-Schwartz inequality to derive

where the last step arises from conditions that $\frac{\sigma}{\lambda_{\min}^{\star}}\sqrt{\frac{d\log d}{p}}\ll\frac{1}{\mu^{3/2}r}$ and $\kappa\asymp 1$.

(5) Taking (126), (129), (132) and (133) together, we can invoke the sample size assumption that $p\gg\mu^{3}r^{2}d^{-2}\log^{3}d$ and the union bound to show that: with probability greater than $1-O(d^{-10})$ one has

provided 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.

A.5 Proof of Lemma 5.5

Fix an arbitrary $m\in[d]$. Recall our notation of $\bm{\Delta}^{t,\left(m\right)}_{\bm{T}}$ in (124). To simplify presentation, we further define

leaving us with two terms to deal with. As it turns out, we will show that $\alpha_{1}$ is the dominant term and $\alpha_{2}$ is negligible. To simplify presentation, we shall assume that $\{E_{i,j,k}\}_{i,j,k\in[d]}$ (resp. $\{\chi_{i,j,k}\}_{i,j,k\in[d]}$) are independent random variables.

Regarding $\alpha_{1}$, the definition of $\mathcal{P}_{\Omega_{-m}}$ allows us to derive

We can express \bm{\Delta}^{t,\left(m\right)}_{\bm{T}}=\sum_{s\in[r]}\big{(}\bm{\Delta}_{s}^{t,\left(m\right)}+\bm{u}_{s}^{\star}\big{)}^{\otimes 3}-\bm{u}_{s}^{\star\otimes 3} and compute that

for each $s\in[r]$. This further indicates that

In what follows, we will control the four terms separately.

For $\beta_{1}$, by (119), we use Cauchy-Schwarz to show that

Regarding $\beta_{2}$, by (119), we apply the Cauchy-Schwarz inequality again to get that: for each $s\in[r]$,

Together with the assumption that $\kappa\asymp 1$, this implies that

With regards to $\beta_{3}$, we show that for each $s\in[r]$:

Turning attention to $\beta_{4}$, we observe that for each $s\in[r]$,

where the last line follows from (119) and the rank assumption $r\ll\sqrt{d/\mu}$. As a consequence,

With regards to $\alpha_{2}$, it follows from the definition (46) and (134) that

Recall the definition of $\mathcal{P}_{\Omega_{-m}}$ and $\mathcal{P}_{m}$. For the $m$-th row, we have

From the triangle inequality, we can further decompose

Let us consider $\gamma_{1}$ first. It is straightforward to calculate that

for each $s\in[r]$. From (119), we use the triangle inequality and the Cauchy-Schwarz inequality to upper bound

Moreover, it is easy to see that the upper bound also holds for $\gamma_{2}$. As for $\gamma_{3}$, we can express

Similarly, we combine (119) wtih the triangle inequality and the Cauchy-Schwarz inequality to bound

where the last inequality follows from the assumption that $\kappa\asymp 1$.

Putting (141) and (145) together, we reach the conclusion from (48) and the condition $r\ll\sqrt{d/\mu}$ that,

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

Recognizing that the above bound holds for any $1\leq m\leq d$, we conclude the proof.

A.6 Proof of Lemma 5.6

Combining Lemma 5.4 and Lemma 5.5, we conclude that with probability at least $1-O(d^{-10})$,

with the proviso that $C_{3}/\left(C_{5}+C_{7}\right)$ and $C_{4}/\left(C_{6}+C_{8}\right)$ are both sufficiently large.

Appendix B Proofs for retrieving tensor components

We shall often operate upon the event where the claims in Lemma 5.7 hold, which happens with very high probability (i.e. at least $1-O\left(d^{-10}\right)$). Recall the definition of $\bm{\gamma}^{\star\tau}$ in (72). Since $\bm{\theta}^{\tau}=\bm{U}\bm{U}^{\top}\bm{g}^{\tau}$, this allows us to write that: for each $1\leq i\leq r$,

where we recall that $\lambda_{i}^{\star}=\|\bm{u}_{i}^{\star}\|_{2}^{3}$ and $\overline{\bm{u}}_{i}^{\star}=\bm{u}_{i}^{\star}/\|\bm{u}_{i}^{\star}\|_{2}$. Given that $\bm{g}^{\tau}$ is a Gaussian vector independent of $\bm{U}$, we observe that $\bm{\gamma}^{\star\tau}$ is zero-mean Gaussian conditional on $\Omega$ and $\bm{E}$. In order to understand the order statistics associated with this vector, we first look at its covariance matrix.

Denote by $\bm{\Sigma}^{\tau}$ the covariance matrix of $\bm{\gamma}^{\star\tau}$ (conditional on $\bm{U}$). Then we have

where we denote by $\mathcal{P}_{\bm{U}}(\bm{z})=\bm{U}\bm{U}^{\top}\bm{z}$ and $\mathcal{P}_{\bm{U}^{\perp}}(\bm{z})=\left(\bm{I}-\bm{U}\bm{U}^{\top}\right)\bm{z}$. In addition, since the unit vector $\overline{\bm{u}}_{i}^{\star}$ lies in the span of the columns of $\bm{U}_{\mathsf{orth}}^{\star}$ (cf. (54)), it follows from Lemma D.6 that

where $\bm{R}$ is a rotation matrix defined in (55). This together with (147a) gives

where we have also used the fact that \left\|\mathcal{P}_{\bm{U}}\big{(}\overline{\bm{u}}_{i}^{\star}\big{)}\right\|_{2}\leq\left\|\overline{\bm{u}}_{i}^{\star}\right\|_{2}=1. Moreover, taking together (147a), (148b) and the incoherence condition, we see that for any $1\leq i\neq j\leq r$,

which is expected to be small if $\delta_{1}$ is small.

From our assumptions on the sample size, the rank and the condition number, we can invoke Lemma 5.7 to see that $\kappa r\delta_{1}\ll 1$, where $\delta_{1}$ is defined in (150) and $\kappa=\lambda_{\max}^{\star}/\lambda_{\min}^{\star}$. Thus, we can decompose $\bm{\Sigma}^{\tau}$ into two components as follows

As it turns out, both $\widehat{\bm{\Sigma}}^{\tau}$ and $\breve{\bm{\Sigma}}^{\tau}$ are positive definite. Indeed, we first learn from (149) and (150) that: the $i$-th digaonal entry of $\breve{\bm{\Sigma}}^{\tau}$ obeys

where (i) holds since $\delta_{1}<1$ under our assumptions, (ii) follows since $\kappa\lambda_{i}^{\star}\geq\kappa\lambda_{\min}^{\star}=\lambda_{\max}^{\star}$, and the last line makes use of (150). This implies that $\breve{\bm{\Sigma}}^{\tau}$ is diagonally dominant, and hence $\breve{\bm{\Sigma}}^{\tau}\succeq\bm{0}$. In conclusion, both $\widehat{\bm{\Sigma}}^{\tau}$ and $\breve{\bm{\Sigma}}^{\tau}$ are positive definite.

Let $\widehat{\bm{\gamma}}^{\star\tau}$ and $\breve{\bm{\gamma}}^{\star\tau}$ be independent zero-mean Gaussian random vectors with covariance matrices $\widehat{\bm{\Sigma}}^{\tau}$ and $\breve{\bm{\Sigma}}^{\tau}$, respectively. Clearly, the distribution of $\bm{\gamma}^{\star\tau}$ is identical to that of $\widehat{\bm{\gamma}}^{\star\tau}+\breve{\bm{\gamma}}^{\star\tau}$. Consequently, it allows us to look at the distributions of these two random vectors separately.

In view of (149) and the fact $\kappa r\delta_{1}<1$, one has

Thus, with probability at least $1-O\left(d^{-10}\right)$, we have

We then turn attention to $\widehat{\bm{\gamma}}^{\star\tau}$, which is composed of independent Gaussian random variables. Let us define \widehat{\Delta}^{\tau}:=\widehat{\gamma}_{1}^{\star\tau}-\max_{1<i\leq r}\big{|}\widehat{\gamma}_{i}^{\star\tau}\big{|} for each $1\leq\tau\leq L$ and let $\widehat{\Delta}_{1}^{\left(1\right)}\geq\widehat{\Delta}_{1}^{\left(2\right)}\geq\dots\geq\widehat{\Delta}_{1}^{\left(L\right)}$ denote the order statistics of \big{\{}\widehat{\Delta}_{1}^{\tau}\big{\}}_{\tau=1}^{L} in descending order. Fix any small constant $\delta>0$. Invoke Lemma D.5 to demonstrate that: with probability greater than $1-\delta/r$,

Putting (152) and (154) together and invoking the triangle inequality immediately establish (74a) and (74b). On the other hand, combining (152) with (155a) proves (75a); (153a) and (155b) taken collectively establish (75b), whereas (153b) and (155c) prove (75c).

B.2 Proof of Lemma 5.13

Recall that the vector of interest $\overline{\bm{u}}^{\tau}$ is the leading singular vector of $\bm{M}^{\tau}$ (as constructed in (70c)), where $\bm{M}^{\tau}$ satisfies

Instate the assumptions of Lemma 5.13. With probability at least $1-O\left(d^{-10}\right)$, one has

As a consequence, recalling the definition of $\mathcal{E}_{\mathsf{op}}$ and $\mathcal{E}_{\mathsf{proj}}$ in (81) and (79) respectively, one has

It then follows from Weyl’s inequality that

where $\sigma_{i}(\bm{Z})$ denotes the $i$-th largest singular value of a matrix $\bm{Z}$ and we use the condition that $\mathcal{E}_{\mathsf{op}}\ll 1$. All in all, these arguments justify that the spectrum of $\bm{M}^{\tau}$ is fairly close to that of $\bm{M}^{\star\tau}$.

Next, we look at the gap between the two leading singular values of $\bm{M}^{\star\tau}$. To begin with, it is self-evident from the definition of $\bm{M}^{\star\tau}$ that: $\overline{\bm{u}}_{1}^{\star}$ is the singular vector of $\bm{M}^{\star\tau}$. In fact, we claim one further result, that is, $\overline{\bm{u}}_{1}^{\star}$ is indeed the leading singular vector of $\bm{M}^{\star\tau}$ whose singular value is given by \sigma_{1}\big{(}\bm{M}^{\star\tau}\big{)}=\gamma_{1}^{\star\tau}. Towards this end, let us define

Let $|\gamma^{\star\tau}|_{(1)}\geq\dots\geq|\gamma^{\star\tau}|_{(r)}$ denote the absolute values of $\{\gamma_{i}^{\star\tau}\}_{i=1}^{r}$ in descending order. This together with Lemma 5.12 implies that

as long as $\kappa r\sqrt{(\mu\log d)/d}\ll 1$. Given that $\overline{\bm{u}}_{1}^{\star}$ is the singular vector of $\bm{M}^{\star\tau}$ with singular value $\gamma_{1}^{\star\tau}$, we can conclude that \sigma_{1}\big{(}\bm{M}^{\star\tau}\big{)}=\gamma_{1}^{\star\tau}. This also allows us to lower bound the gap between the two largest singular values $\bm{M}^{\star\tau}$ as follows

provided that $\kappa r\sqrt{(\mu\log d)/d}\ll 1$. We also know from (163) and (164) that

Combined with (162) and Wedin’s theorem, we conclude that

Here, we have made use of the fact that $\overline{\bm{u}}^{\tau}$ is the leading singular vector of $\bm{M}^{\tau}$ obeying \big{\langle}\overline{\bm{u}}^{\tau},\overline{\bm{u}}_{1}^{\star}\big{\rangle}\geq 0.

B.3 Proof of Lemma B.1

We first consider the spectral norm of $\bm{F}^{\tau}$. Recall the definition that $\bm{\theta}^{\tau}=\bm{U}\bm{U}^{\top}\bm{g}^{\tau}$. Let us define

In the sequel, we shall control these two terms separately.

To bound $\left\|\bm{X}\right\|$, observe that $\bm{\theta}^{\star\tau}$ is independent of $p^{-1}\bm{T}-\bm{T}^{\star}$. By Lemma D.4, one has

Combining (166) with (167) and (168) reveals that with probability exceeding $1-O\left(d^{-20}\right)$,