Nonconvex Matrix Factorization from Rank-One Measurements

Introduction

Quantum state tomography. Estimating the density operator of a quantum system can be formulated as a low-rank positive semidefinite matrix recovery problem using rank-one measurements, when the density operator is almost pure . A problem of similar mathematical formulation occurs in phase space tomography , where the goal is to reconstruct the correlation function of a wave field.

Due to the quadratic nature of the measurements, the natural least-squares empirical risk formulation is highly nonconvex and in general challenging to solve. To be more specific, consider the following optimization problem:

which aims to optimize a degree-4 polynomial in ${\bm{X}}$ and is NP hard in general. The problem, however, may become tractable under certain random designs, and may even be solvable using simple methods like gradient descent. Our main finding is the following: under i.i.d. Gaussian design (i.e. ${\bm{a}}_{i}\sim\mathcal{N}(\bm{0},{\bm{I}}_{n})$), vanilla gradient descent combined with spectral initialization achieves appealing performance guarantees both statistically and computationally.

Statistically, we show that gradient descent converges exactly to the true factor ${\bm{X}}^{\natural}$ (modulo unrecoverable global ambiguity), as soon as the number of measurements exceeds the order of $O(nr^{4}\log n)$. When $r$ is fixed independent of $n$, this sample complexity is near-optimal up to some logarithmic factor.

Computationally, to achieve $\epsilon$-accuracy, gradient descent requires an iteration complexity of $O(r^{2}\log(1/\epsilon))$ (up to logarithmic factors), with a per-iteration cost of $O(mnr)$. When $r$ is fixed independent of $m$ and $n$, the computational complexity scales linearly with $mn$, which is proportional to the time taken to read all data.

These findings significantly improve upon existing results that require either resampling (which is not sample-efficient and is not the algorithm one actually runs in practice ), or high iteration complexity (which results in high computation cost ). In particular, our work is most related to that also studied the effectiveness of gradient descent. The results in require a sample complexity on the order of $nr^{6}\log^{2}{n}$, as well as an iteration complexity of $O(n^{4}r^{2}\log(1/\epsilon))$ (up to logarithmic factors) to attain $\epsilon$-accuracy. In comparison, our theory improves the sample complexity to $O(nr^{4}\log n)$ and, perhaps more importantly, establishes a much lower iteration complexity of $O(r^{2}\log(1/\epsilon))$ (up to logarithmic factor). To the best of our knowledge, this work is the first nonconvex algorithm (without resampling) that achieves both near-optimal statistical and computational guarantees with respect to $n$.

2 Surprising Effectiveness of Gradient Descent

Recently, gradient descent has been widely employed to address various nonconvex optimization problems due to its appealing efficiency from both statistical and computational perspectives. Despite the nonconvexity of (1), showed that within a local neighborhood of ${\bm{X}}^{\natural}$, where ${\bm{X}}$ satisfies

$f({\bm{X}})$ behaves like a strongly convex function, at least along certain descending directions. However, this region itself is not enough to guarantee computational efficiency, and consequently, the smoothness parameter derived in is as large as $n^{2}$ (even ignoring additional polynomial factors in $r$), leading to a step size as small as $O(1/n^{4})$ and an iteration complexity of $O(n^{4}\log(1/\epsilon))$. These are fairly pessimistic.

In order to improve computational guarantees, it might be tempting to employ appropriately designed regularization operations — such as truncation and projection . These explicit regularization operations are capable of stabilizing the search direction, and make sure the whole trajectory is in a basin of attraction with benign curvatures surrounding the ground truth. However, such explicit regularizations complicate algorithm implementations, as they introduce more tuning parameters.

Our work is inspired by , which uncovers the “implicit regularization” phenomenon of vanilla gradient descent for nonconvex estimation problems such as phase retrieval and low-rank matrix completion. In words, even without extra regularization operations, vanilla gradient descent always follows a path within some region around the global optimum with nice geometric structure, at least along certain directions. The current paper demonstrates that a similar phenomenon persists in low-rank matrix factorization from rank-one measurements.

To describe this phenomenon in a precise manner, we need to specify which region enjoys the desired geometric properties. To this end, consider a local region around ${\bm{X}}^{\natural}$ where ${\bm{X}}$ is “incoherent”This is called incoherent because if ${\bm{X}}$ is aligned (and hence coherent) with the sensing vectors, \big{\|}{\bm{a}}_{l}^{\top}\big{(}{\bm{X}}-{\bm{X}}^{\natural}\big{)}\big{\|}_{2} can be $O(\sqrt{n})$ times larger than the right-hand side of (3). with all sensing vectors in the following sense:

We term the intersection of (2) and (3) the Region of Incoherence and Contraction (RIC). The nice feature of the RIC is this: within this region, the loss function $f({\bm{X}})$ enjoys a smoothness parameter that scales as $O(\max\{r,\log n\})$ (namely, $\|\nabla^{2}f({\bm{x}})\|\lesssim\max\{r,\log n\}$, which is much smaller than $O(n^{2})$ provided in ). As is well known, a region enjoying a smaller smoothness parameter enables more aggressive progression of gradient descent.

A key question remains as to how to prove that the trajectory of gradient descent never leaves the RIC. This is, unfortunately, not guaranteed by standard optimization theory, which only ensures contraction of the Euclidean error. To address this issue, we resort to the leave-one-out trick that produces auxiliary trajectories of gradient descent that use all but one sample. This allows us to establish the incoherence condition by leveraging the statistical independence of the leave-one-out trajectory w.r.t. the corresponding sensing vector that has been left out. Our theory refines the leave-one-out argument and further establishes linear contraction in terms of the entry-wise prediction error.

3 Notations

Algorithms and Main Results

To begin with, we present the formal problem setup. Suppose we are given a set of $m$ rank-one measurements as follows

Throughout this paper, we assume the condition number is bounded by some constant independent of $n$ and $r$, i.e. $\kappa=O(1)$. Our goal is to recover ${\bm{X}}^{\natural}$, up to (unrecoverable) orthonormal transformation, from the measurements ${\bm{y}}=\left\{y_{i}\right\}_{i=1}^{m}$ in a statistically and computationally efficient manner.

The algorithm studied herein is a combination of vanilla gradient descent and a judiciously designed spectral initialization. Specifically, consider minimizing the squared loss:

which is a nonconvex function. We attempt to optimize this function iteratively via gradient descent

where ${\bm{X}}_{t}$ denotes the estimate in the $t$th iteration, $\mu_{t}$ is the step size/learning rate, and the gradient $\nabla f({\bm{X}})$ is given by

For initialization, similar to ,Compared with , when setting the eigenvalues in (10), we use the sample mean $\lambda$ rather than $\lambda_{r+1}\left({\bm{Y}}\right)$ to estimate $\frac{1}{2}\|{\bm{X}}^{\natural}\|_{{\mathsf{F}}}^{2}$. we apply the spectral method, which sets the columns of ${\bm{X}}_{0}$ as the top-$r$ eigenvectors — properly scaled — of a matrix ${\bm{Y}}$ as defined in (9). The rationale is this: the mean of ${\bm{Y}}$ is given by

and hence the principal components of ${\bm{Y}}$ form a reasonable estimate of ${\bm{X}}^{\natural}$, provided that there are sufficiently many samples. The full algorithm is described in Algorithm 1.

2 Performance Guarantees

with $\mathcal{O}^{r\times r}$ denoting the set of all $r\times r$ orthonormal matrices. Accordingly, we have the following theoretical performance guarantees of Algorithm 1.

Suppose that $m\geq cnr^{3}(r+\sqrt{\kappa})\kappa^{3}\log{n}$ with some large enough constant $c>0$, and that the step size obeys $0<\mu_{t}:=\mu=\frac{c_{4}}{(r\kappa+\log{n})^{2}\sigma_{r}^{2}({\bm{X}}^{\natural})}$. Then with probability at least $1-O(mn^{-7})$, the iterates satisfy

for all $0\leq t\leq c_{3}n^{5}$. Here, $c_{1},\cdots,c_{4}$ are some universal positive constants.

The precise expression of required sample complexity in Theorem 1 can be written as $m\geq c\max\left\{\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}{\sigma_{r}\left({\bm{X}}^{\natural}\right)}\sqrt{r},\kappa\right\}\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{5}}{\sigma_{r}^{5}\left({\bm{X}}^{\natural}\right)}n\sqrt{r}\log{\left(\kappa n\right)}$ with some large enough constant $c>0$. By adjusting constants, with probability at least $1-O(mn^{-7})$, (15) holds for $0\leq t\leq O(n^{c_{5}})$ in any power $c_{5}\geq 1$.

Theorem 1 has the following implications.

Near-optimal sample complexity when $r$ is fixed: Theorem 1 suggests that spectrally-initialized vanilla gradient descent succeeds as soon as $m=O(nr^{4}\log n)$. When $r=O(1)$, this leads to near-optimal sample complexity up to logarithmic factor. In fact, once the spectral initialization is finished, a sample complexity at $m=O(nr^{3}\log n)$ can guarantee the linear convergence to the global optima. To the best of our knowledge, this outperforms all performance guarantees in the literature obtained for any nonconvex method without requiring resampling.

Implicit regularization: Theorem 1 demonstrates that both the spectral initialization and the gradient descent updates provably control the entry-wise error \max_{1\leq l\leq m}\left\|{\bm{a}}_{l}^{\top}\big{(}{\bm{X}}_{t}{\bm{Q}}_{t}-{\bm{X}}^{\natural}\big{)}\right\|_{2}, and the iterates remain incoherent with respect to all the sensing vectors. In fact, the entry-wise error decreases linearly as well, which is not characterized in .

Theorem 1 is established using a fixed step size. According to our theoretical analysis, the incoherence condition (15) has a significant impact on the convergence rate. After a few iterations, the incoherence condition can be bounded independent of $\log n$, which leads to a larger step size and faster convergence. Specifically, we have the following corollary.

Under the same setting of Theorem 1, after $T_{a}=c_{6}\max\{\kappa^{2}r^{2}\log{n},\log^{3}{n}\}$ iterations, the step size can be relaxed as $0<\mu_{t}:=\mu=\frac{c_{7}}{r^{2}\kappa^{2}\sigma_{r}^{2}({\bm{X}}^{\natural})}$, with some universal constant $c_{6},c_{7}>0$, then the iterates satisfy

for all $t\geq T_{a}$, with probability at least $1-O(mn^{-7})$.

Related Work

Instead of directly estimating ${\bm{X}}^{\natural}$, the problem of interest can be also solved by estimating ${\bm{M}}^{\natural}={\bm{X}}^{\natural}{\bm{X}}^{\natural\top}$ in higher dimension via nuclear norm minimization, which requires $O(nr)$ measurements for exact recovery . See also for the phase retrieval problem. However, nuclear norm minimization, often cast as the semidefinite programming, is in general computationally expensive to deal with large-scale data.

On the other hand, nonconvex approaches have drawn intense attention in the past decade due to their ability to achieve computational and statistical efficiency all at once. Specifically, for the phase retrieval problem, Wirtinger Flow (WF) and its variants have been proposed. As a two-stage algorithm, it consists of spectral initialization and iterative gradient updates. This strategy has found enormous success in solving other problems such as low-rank matrix recovery and completion , blind deconvolution , and spectral compressed sensing . We follow a similar route but analyze a more general problem that includes phase retrieval as a special case.

The paper is most close to our work, which studied the local convexity of the same loss function and developed performance guarantees for gradient descent using a similar, but different spectral initialization scheme. As discussed earlier, due to the pessimistic estimate of the smoothness parameter, they only allow a diminishing learning rate (or step size) of $O(1/n^{4})$, leading to a high iteration complexity. We not only provide stronger computational guarantees, but also improve the sample complexity, compared with .

Several other existing works have suggested different approaches for low-rank matrix factorization from rank-one measurements, of which the statistical and computational guarantees to reach $\epsilon$-accuracy are summarized in Table 1. We note our guarantee is the only one that achieves simultaneous near-optimal sample complexity and computational complexity. Iterative algorithms based on alternating minimization or noisy power iterations require a fresh set of samples at every iteration, which is never executed in practice, and the sample complexity grows unbounded for exact recovery.

Many nonconvex methods have been proposed and analyzed recently to solve the phase retrieval problem, including the Kaczmarz method and approximate message passing . In , the Kaczmarz method is generalized to solve the problem studied in this paper, but no theoretical performance guarantees are provided.

The local geometry studied in our paper is in contrast to , which studied the global landscape of phase retrieval, and showed that there are no spurious local minima as soon as the sample complexity is above $O(n\log^{3}n)$. It will be interesting to study the landscape property of the generalized model in our paper.

Our model is also related to learning shallow neural networks. studied the performance of gradient descent with resampling and an initialization provided by the tensor method for various activation functions, however their analysis did not cover quadratic activations. For quadratic activations, adopts a greedy learning strategy, and can only guarantee sublinear convergence rate. Moreover, studied the optimization landscape for an over-parameterized shallow neural network with quadratic activation, where $r$ is larger than $n$.

Outline of Theoretical Analysis

This section provides the proof sketch of the main results, with the details deferred to the appendix. Our theoretical analysis is inspired by the work of for phase retrieval and follows the general recipe outlined in , while significant changes and elaborate derivations are needed. We refine the analysis to show that both the signal reconstruction error and the entry-wise error contract linearly, where the latter is not revealed by . In below, we first characterize a region of incoherence and contraction that enjoys both strong convexity and smoothness along certain directions. We then demonstrate — via an induction argument — that the iterates always stay within this nice region. Finally, the proof is complete by validating the desired properties of spectral initialization.

We start with characterizing a local region around ${\bm{X}}^{\natural}$, within which the loss function enjoys desired restricted strong convexity and smoothness properties. This requires exploring the property of the Hessian of $f({\bm{X}})$, which is given by

Suppose the sample size obeys $m\geq c\frac{\|{\bm{X}}^{\natural}\|_{{\mathsf{F}}}^{4}}{\sigma_{r}^{4}\left({\bm{X}}^{\natural}\right)}nr\log{\left(n\kappa\right)}$ for some sufficiently large constant $c>0$. Then with probability at least $1-c_{1}n^{-12}-me^{-1.5n}-mn^{-12}$, we have

hold simultaneously for all matrices ${\bm{X}}$ and ${\bm{V}}$ satisfying the following constraints:

and ${\bm{V}}={\bm{T}}_{1}{\bm{Q}}_{{\bm{T}}}-{\bm{T}}_{2}$ satisfying

The condition (20) on ${\bm{X}}$ formally characterizes the RIC, which enjoys the claimed restricted strong convexity (see (18)) and smoothness (see (19)). With Lemma 1 in mind, it is easy to see that if ${\bm{X}}_{t}$ lies within the RIC, the estimation error shrinks in the presence of a properly chosen step size. This is given in the lemma below whose proof can be found in Appendix D.

Suppose the sample size obeys $m\geq c\frac{\|{\bm{X}}^{\natural}\|_{{\mathsf{F}}}^{4}}{\sigma_{r}^{4}\left({\bm{X}}^{\natural}\right)}nr\log{\left(n\kappa\right)}$ for some sufficiently large constant $c>0$. Then with probability at least $1-c_{1}n^{-12}-me^{-1.5n}-mn^{-12}$, if ${\bm{X}}_{t}$ falls within the RIC as described in (20), we have

provided that the step size obeys 0<\mu_{t}\equiv\mu\leq\frac{1.026\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\big{(}1.5\sigma_{r}^{2}({\bm{X}}^{\natural})\log{n}+6\|{\bm{X}}^{\natural}\|_{{\mathsf{F}}}^{2}\big{)}^{2}}. Here, $c_{1}>0$ is some universal constant.

Assuming that the iterates $\{{\bm{X}}_{t}\}$, stay within the RIC (see (20)) for the first $T_{c}$ iterations, according to Lemma 2, we have, by induction, that

for some large enough constant $c$. The iterates when $t\geq T_{c}$ are easier to deal with; in fact, it is easily seen that ${\bm{X}}_{t+1}$ stays in the RIC since

2 Introducing Leave-One-Out Sequences

It has now become clear that the key remaining step is to verify that the iterates $\{{\bm{X}}_{t}\}$ satisfy (20) for the first $T_{c}$ iterations, where $T_{c}$ is on the order of (22). Verifying (20b) is conceptually hard since the iterates $\{{\bm{X}}_{t}\}$ are statistically dependent with all the sensing vectors $\{{\bm{a}}_{i}\}_{i=1}^{m}$. To tackle this problem, for each $1\leq l\leq m$, we introduce an auxiliary leave-one-out sequence $\{{\bm{X}}^{(l)}_{t}\}$, which discards a single measurement from consideration. Specifically, the sequence $\{{\bm{X}}^{(l)}_{t}\}$ is the gradient iterates operating on the following leave-one-out function

See Algorithm 2 for a formal definition of the leave-one-out sequences. Again, we want to emphasize that Algorithm 2 is just an auxiliary procedure useful for the theoretical analysis, and it does not need to be implemented in practice.

3 Establishing Incoherence via Induction

Our proof is inductive in nature with the following induction hypotheses:

Furthermore, the step size $\mu$ is chosen as

with appropriate universal constant $c_{0}>0$.

Our goal is to show that if the $t$th iteration ${\bm{X}}_{t}$ satisfies the induction hypotheses (28), then the $(t+1)$th iteration ${\bm{X}}_{t+1}$ also satisfies (28). It is straightforward to see that the hypothesis (28a) has already been established by Lemma 2, and we are left with (28b) and (28c). We first establish (28b) in the following lemma, which measures the proximity between ${\bm{X}}_{t}$ and the leave-one-out versions ${\bm{X}}_{t}^{(l)}$, whose proof is provided in Appendix E.

Suppose the sample size obeys $m\geq c\frac{\|{\bm{X}}^{\natural}\|_{{\mathsf{F}}}^{4}}{\sigma_{r}^{4}\left({\bm{X}}^{\natural}\right)}nr\log{\left(n\kappa\right)}$ for some sufficiently large constant $c>0$. If the induction hypotheses (28) hold for the $t$th iteration, with probability at least $1-c_{1}n^{-12}-me^{-1.5n}-mn^{-12}$, we have

as long as the step size obeys (30). Here, $c_{1}>0$ is some absolute constant.

In addition, the incoherence property of ${\bm{X}}_{t+1}^{(l)}$ with respect to the $l$th sensing vector ${\bm{a}}_{l}$ is relatively easier to establish, due to their statistical independence. Combined with the proximity bound from Lemma 3, this allows us to justify the incoherence property of the original iterates ${\bm{X}}_{t+1}$, as summarized in the lemma below, whose proof is given in Appendix F.

Suppose the sample size obeys $m\geq c\frac{\|{\bm{X}}^{\natural}\|_{{\mathsf{F}}}^{4}}{\sigma_{r}^{4}\left({\bm{X}}^{\natural}\right)}nr\log{\left(n\kappa\right)}$ for some sufficiently large constant $c>0$. If the induction hypotheses (28) hold for the $t$th iteration, with probability exceeding $1-c_{1}n^{-12}-me^{-1.5n}-2mn^{-12}$,

holds as long as the step size satisfies (30). Here, $c_{1}>0$ is some universal constant.

4 Spectral Initialization

Finally, it remains to verify that the induction hypotheses hold for the initialization, i.e. the base case when $t=0$. This is supplied by the following lemma, whose proof is given in Appendix G.

Suppose that the sample size exceeds $m\geq c\max\left\{\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}{\sigma_{r}\left({\bm{X}}^{\natural}\right)}\sqrt{r},\kappa\right\}\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{5}}{\sigma_{r}^{5}\left({\bm{X}}^{\natural}\right)}n\sqrt{r}\log{n}$ for some sufficiently large constant $c>0$. Then ${\bm{X}}_{0}$ satisfies (28) with probability at least $1-c_{1}n^{-12}-me^{-1.5n}-3mn^{-12}$, where $c_{1}$ is some absolute positive constant.

Conclusions

In this paper, we have shown that low-rank positive semidefinite matrices can be recovered from a near-minimal number of random rank-one measurements, via the vanilla gradient descent algorithm following spectral initialization. Our results significantly improve upon existing results in several ways, both computationally and statistically. In particular, our algorithm does not require resampling at every iteration (and hence requires fewer samples). The gradient iteration can provably employ a much more aggressive step size than what was suggested in prior literature (e.g. ), thus resulting in much smaller iteration complexity and hence lower computational cost. All of this is enabled by establishing the implicit regularization feature of gradient descent for nonconvex statistical estimation, where the iterates remain incoherent with the sensing vectors throughout the execution of the whole algorithm.

There are several problems that are worth exploring in future investigation. For example, our theory reveals the typical size of the fitting error of ${\bm{X}}_{t}$ (i.e. $y_{i}-\|{\bm{a}}_{i}^{\top}{\bm{X}}^{\natural}\|_{2}$) in the presence of noiseless data, which would serve as a helpful benchmark when separating sparse outliers in the more realistic scenario. Another direction is to explore whether implicit regularization remains valid for learning shallow neural networks . Since the current work can be viewed as learning a one-hidden-layer fully-connected network with a quadratic activation function $\sigma(z)=z^{2}$, it would be of great interest to study if the techniques utilized herein can be used to develop strong guarantees when the activation function takes other forms.

Acknowledgements

The work of Y. Li and Y. Chi is supported in part by AFOSR under the grant FA9550-15-1-0205, by ONR under the grant N00014-18-1-2142, and by NSF under the grants CAREER ECCS-1818571 and CCF-1704245.

Appendices

Appendix A Technical Lemmas

In this section, we document a few useful lemmas that are used throughout the proof.

where $\Phi(\cdot)$ is the cumulative distribution function of a standard Gaussian variable.

[39, Theorem 5.39] Suppose the ${\bm{a}}_{i}$’s are i.i.d. random vectors following ${\bm{a}}_{i}\sim\mathcal{N}\left(\bm{0},{\bm{I}}_{n}\right)$, $i=1,\cdots,m$. Then for every $t\geq 0$ and $0<\delta\leq 1$,

Suppose the ${\bm{a}}_{i}$’s are i.i.d. random vectors following ${\bm{a}}_{i}\sim\mathcal{N}\left(\bm{0},{\bm{I}}_{n}\right)$, $i=1,\cdots,m$. Then with probability at least $1-me^{-1.5n}$, we have

Define ${\bm{W}}=[{\bm{w}}_{1},{\bm{w}}_{2},\cdots,{\bm{w}}_{r}]$, then we can write $\left\|{\bm{a}}_{i}^{\top}{\bm{W}}\right\|_{2}^{2}=\sum_{k=1}^{r}\left({\bm{a}}_{i}^{\top}{\bm{w}}_{k}\right)^{2}$. Recognize that $\left({\bm{a}}_{i}^{\top}\frac{{\bm{w}}_{k}}{\left\|{\bm{w}}_{k}\right\|_{2}}\right)^{2}$ follows the $\chi^{2}$ distribution with $1$ degree of freedom. It then follows from [40, Lemma 1] that

for any $t>0$. Taking $t=13\log{n}$ yields

Finally, taking the union bound, we obtain

with probability at least $1-mrn^{-13}$. ∎

Let ${\bm{X}}=[{\bm{x}}_{1},{\bm{x}}_{2},\cdots,{\bm{x}}_{r}]$ and ${\bm{H}}=[{\bm{h}}_{1},{\bm{h}}_{2},\cdots,{\bm{h}}_{r}]$. Based on the simple facts

with probability at least $1-c_{1}rn^{-13}$, where $c_{1}>0$ is some absolute constant.

This proof adapts the results of [2, Lemma 7.4] with refining the probabilities. Let ${\bm{a}}(1)$ be the first element of a vector ${\bm{a}}\sim\mathcal{N}\left(\bm{0},{\bm{I}}_{n}\right)$. Based on [41, Theorem 1.9], we have

So, by setting $m\gg\delta^{-2}n$, we have

with probability at least $1-c_{4}n^{-13}$ for some constant $c_{4}>0$. Moreover, following [40, Lemma 1], we know

if setting $t=14\log{m}$. Therefore, as long as $m\geq cn$, we have

with probability at least $1-c_{5}n^{-13}$ for some constant $c_{5}>0$.

With (31) and (32), the results in [2, Lemma 7.4] imply that for any $0<\delta\leq 1$, as soon as $m\geq c\delta^{-2}n\log{n}$ for some sufficiently large constant $c$, with probability at least $1-c_{1}n^{-13}$,

Appendix B Proof of Lemma 1

The crucial ingredient for proving the lower bound (18) is the following lemma, whose proof is provided in Appendix C.

Suppose m\geq c\frac{\big{\|}{\bm{X}}^{\natural}\big{\|}_{{\mathsf{F}}}^{4}}{\sigma_{r}^{4}\left({\bm{X}}^{\natural}\right)}nr\log{\left(n\kappa\right)} with some large enough positive constant $c$, then with probability at least $1-c_{1}n^{-12}-me^{-1.5n}$, we have

for all matrices ${\bm{X}}$ and ${\bm{V}}$ where ${\bm{X}}$ satisfies $\left\|{\bm{X}}-{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}\leq\frac{1}{24}\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}$. Here, $c_{1}>0$ is some universal constant.

With Lemma 14 in place, we are ready to prove (18). Let ${\bm{V}}={\bm{T}}_{1}{\bm{Q}}_{{\bm{T}}}-{\bm{T}}_{2}$ satisfy the assumptions in Lemma 1, then we can demonstrate that

To prove the upper bound (19) asserted in the lemma, we make the observation that the Hessian in (17) satisfies

where (37) follows from the fact $\left\|{\bm{I}}\otimes{\bm{A}}\right\|=\left\|{\bm{A}}\right\|$. It is seen from Lemma 13 that

when setting \delta\leq 0.02\frac{\sigma_{r}^{2}\big{(}{\bm{X}}^{\natural}\big{)}}{\big{\|}{\bm{X}}^{\natural}\big{\|}_{{\mathsf{F}}}^{2}}. Moreover, it is straightforward to check that

With regards to the first term $B_{1}$, note that by Lemma 11 and (20b), we can bound

where the last line follows from Lemma 9. The proof is then finished by combining (38) with the preceding bounds on $B_{1}$, $B_{2}$ and $B_{3}$.

Appendix C Proof of Lemma 14

Without loss of generality, we assume $\left\|{\bm{V}}\right\|_{{\mathsf{F}}}=1$. Write

In what follows, we let ${\bm{X}}={\bm{X}}^{\natural}+t\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}{\bm{H}}$ with $t\leq 1/24$ and $\left\|{\bm{H}}\right\|_{{\mathsf{F}}}=1$ which immediately obeys $\left\|{\bm{X}}-{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}\leq\frac{1}{24}\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}$, and express the right-hand side of (40) as

The aim is thus to control $p\left({\bm{V}},{\bm{H}},t\right)$ for all matrices satisfying $\left\|{\bm{H}}\right\|_{{\mathsf{F}}}=1$ and $\left\|{\bm{V}}\right\|_{{\mathsf{F}}}=1$, and for all $t$ obeying $t\leq{1}/{24}$.

We first bound the second term in (41). Let ${\bm{V}}=[{\bm{v}}_{1},{\bm{v}}_{2},\cdots,{\bm{v}}_{r}]$, then by Lemma 13,

By setting $\delta\leq\frac{1}{24}\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{2}}$, we see that with probability at least $1-c_{1}rn^{-13}$,

holds simultaneously for all matrices ${\bm{V}}$, as long as $m\gtrsim\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{4}}{\sigma_{r}^{4}\left({\bm{X}}^{\natural}\right)}n\log{n}$.

Next, we turn to the first term $q\left({\bm{V}},{\bm{H}},t\right)$ in (41), and we need to accommodate all matrices satisfying $\left\|{\bm{H}}\right\|_{{\mathsf{F}}}=1$ and $\left\|{\bm{V}}\right\|_{{\mathsf{F}}}=1$, and all scalars obeying $t\leq{1}/{24}$. The strategy is that we first establish the bound of $q\left({\bm{V}},{\bm{H}},t\right)$ for any fixed ${\bm{H}}$, ${\bm{V}}$ and $t$, and then extend the result to a uniform bound for all ${\bm{H}}$, ${\bm{V}}$ and $t$ by covering arguments.

We will start by assuming that ${\bm{X}}$ and ${\bm{V}}$ are both fixed and statistically independent of $\{{\bm{a}}_{i}\}_{i=1}^{m}$. In view of Lemma 12,

where (45) follows from the Cauchy-Schwarz inequality, (46) comes from the Hölder’s inequality, and (47) is a consequence of Lemma 12. Apply Lemma 8 to arrive at

Substituting ${\bm{X}}={\bm{X}}^{\natural}+t\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}{\bm{H}}$ for ${\bm{X}}$, and using the facts $\left\|{\bm{H}}\right\|_{{\mathsf{F}}}=1$, $\left\|{\bm{V}}\right\|_{{\mathsf{F}}}=1$ and $t\leq 1/24$, we can calculate the following bounds:

C.2 Covering Arguments

Since we have obtained a lower bound on $q\left({\bm{V}},{\bm{H}},t\right)$ for fixed ${\bm{V}}$, ${\bm{H}}$ and $t$, we now move on to extending it to a uniform bound that covers all ${\bm{V}}$, ${\bm{H}}$ and $t$ simultaneously. Towards this, we will invoke the $\epsilon$-net covering arguments for all ${\bm{V}}$, ${\bm{H}}$ and $t$, respectively, and will rely on the fact $\max_{1\leq i\leq m}\left\|{\bm{a}}_{i}\right\|_{2}\leq\sqrt{6n}$ asserted in Lemma 10. For notational convenience, we define

First, consider the $\epsilon$-net covering argument for ${\bm{V}}$. Suppose ${\bm{V}}_{1}$ and ${\bm{V}}_{2}$ are such that $\left\|{\bm{V}}_{1}\right\|_{{\mathsf{F}}}=1$, $\left\|{\bm{V}}_{2}\right\|_{{\mathsf{F}}}=1$, and $\left\|{\bm{V}}_{1}-{\bm{V}}_{2}\right\|_{{\mathsf{F}}}\leq\epsilon$. Then, since

as long as $\epsilon=\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{10584n^{2}\left\|{\bm{X}}^{\natural}\right\|^{2}}$. Based on Lemma 7, the cardinality of this $\epsilon$-net will be

Secondly, consider the $\epsilon$-net covering argument for ${\bm{H}}$. Suppose ${\bm{H}}_{1}$ and ${\bm{H}}_{2}$ obey $\left\|{\bm{H}}_{1}\right\|_{{\mathsf{F}}}=1$, $\left\|{\bm{H}}_{2}\right\|_{{\mathsf{F}}}=1$, and $\left\|{\bm{H}}_{1}-{\bm{H}}_{2}\right\|_{{\mathsf{F}}}\leq\epsilon$. Then one has

as long as $\epsilon=\frac{1}{225n^{2}}\cdot\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}{\left\|{\bm{X}}^{\natural}\right\|}$. Based on Lemma 7, the cardinality of this $\epsilon$-net will be

Finally, consider the $\epsilon$-net covering argument for all $t$, such that $t\leq{1}/{24}$. Suppose $t_{1}$ and $t_{2}$ satisfy $t_{1}\leq{1}/{24}$, $t_{2}\leq{1}/{24}$ and $\left|t_{1}-t_{2}\right|\leq\epsilon$. Then we get

as long as $\epsilon=\frac{1}{5400n^{2}}\cdot\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}{\left\|{\bm{X}}^{\natural}\right\|}$. The cardinality of this $\epsilon$-net will be $\frac{1/24}{\epsilon}\leq cn^{2}\cdot\frac{\left\|{\bm{X}}^{\natural}\right\|}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}$.

Therefore, when m\geq c\frac{\big{\|}{\bm{X}}^{\natural}\big{\|}_{{\mathsf{F}}}^{4}}{\sigma_{r}^{4}\left({\bm{X}}^{\natural}\right)}nr\log{\left(n\kappa\right)} with some large enough constant $c$, for all matrices ${\bm{V}}$ and ${\bm{X}}$ such that $\left\|{\bm{X}}-{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}\leq\frac{1}{24}\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}$, we have

with probability at least $1-e^{-c_{1}nr\log{\left(n\kappa\right)}}-me^{-1.5n}$.

C.3 Finishing the Proof

Appendix D Proof of Lemma 2

Here, (51) follows from the definition of ${\bm{Q}}_{t+1}$ (see (13)), (52) holds owing to the identity $\nabla f\left({\bm{X}}_{t}\right){\bm{Q}}_{t}=\nabla f\left({\bm{X}}_{t}{\bm{Q}}_{t}\right)$ for ${\bm{Q}}_{t}\in\mathcal{O}^{r\times r}$, and (53) arises from the fact that $\nabla f\left({\bm{X}}^{\natural}\right)=\bm{0}$. Let

where $\tau\in$. Then, by the fundamental theorem of calculus for vector-valued functions ,

It is easy to verify that ${\bm{X}}_{t}(\tau)$ satisfies (20) for any $\tau\in$, since

Substituting the above two inequalities into (53) and (55) gives

with the proviso that $\mu\leq\frac{1.026\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left(1.5\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)\log{n}+6\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{2}\right)^{2}}$. This allows us to conclude that

Appendix E Proof of Lemma 3

we will focus on bounding \big{\|}{\bm{X}}_{t+1}{\bm{Q}}_{t}-{\bm{X}}_{t+1}^{(l)}{\bm{R}}_{t}^{(l)}\big{\|}_{{\mathsf{F}}}. Since

we aim to control \big{\|}{\bm{S}}_{t,1}^{(l)}\big{\|}_{{\mathsf{F}}} and \big{\|}{\bm{S}}_{t,2}^{(l)}\big{\|}_{{\mathsf{F}}} separately.

We first bound the term \big{\|}{\bm{S}}_{t,2}^{(l)}\big{\|}_{{\mathsf{F}}}, which is easier to handle. Observe that by Cauchy-Schwarz,

where we have used the triangle inequality, Lemma 10, as well as the induction hypotheses (28c) and (28b). Similarly, the second term in (56) can be bounded as

where we have used (57), Lemma 11, and \sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)\leq\big{\|}{\bm{X}}^{\natural}\big{\|}_{{\mathsf{F}}}^{2}. Similarly, we can also obtain

Substituting (57) and (58) into (56), and using the above inequality, we get

where $C_{4}:=\sqrt{6}C_{3}+C_{2}+11.72$.

Next, we turn to $\left\|{\bm{S}}_{t,1}^{(l)}\right\|_{{\mathsf{F}}}$. By defining

Here, the second line follows from the fundamental theorem of calculus for vector-valued functions , where

for $\tau\in$. Using very similar algebra as in Appendix D, we obtain

It is easy to verify that for all $\tau\in$,

where (62) follows from the induction hypotheses (28a) and (28b), and (63) follows as long as $C_{1}+C_{3}\leq\frac{1}{24}$. Further, for all $1\leq l\leq m$, by the induction hypothesis (28b) and (28c),

as long as $\sqrt{6}C_{3}+C_{2}\leq\frac{1}{24}$. Therefore, Lemma 1 holds for ${\bm{X}}_{t}^{(l)}(\tau)$, and similar to Appendix D, (61) can be further bounded by

as long as $\mu\leq\frac{1.026\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left(1.5\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)\log{n}+6\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{2}\right)^{2}}$. Consequently, combining (59) and (64), we can get

where (65) follows from the induction hypothesis (28b), as long as $m\geq c\kappa\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{2}}{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}n\log{n}$ for some large enough constant $c>0$.

Appendix F Proof of Lemma 4

For any $1\leq l\leq m$, by the statistical independence of ${\bm{a}}_{l}$ and ${\bm{X}}_{t+1}^{(l)}$ and by Lemma 11, we have

as long as $C_{1}\leq\frac{1}{2}$, and following Lemma 3,

as long as $C_{3}\leq\frac{1}{4}$, we can invoke Lemma 37 in and get

Further, by the triangle inequality, Lemma 10, Lemma 3 and Lemma 2, we can deduce that

where the last line follows as long as $5\sqrt{6}C_{3}+5.86C_{1}+29.3C_{3}\leq C_{2}$. The proof is then finished by applying the union bound for all $1\leq l\leq m$.

Appendix G Proof of Lemma 5

then by definition we have ${\bm{Y}}{\bm{Z}}_{0}={\bm{Z}}_{0}\bm{\Sigma}_{0}$, ${\bm{Y}}^{(l)}{\bm{Z}}_{0}^{(l)}={\bm{Z}}_{0}^{(l)}\bm{\Sigma}_{0}^{(l)}$, and

Moreover, let ${\bm{Z}}_{0,c}$ and ${\bm{Z}}_{0,c}^{(l)}$ be the complement matrices of ${\bm{Z}}_{0}$ and ${\bm{Z}}_{0}^{(l)}$, respectively, such that both $\left[{\bm{Z}}_{0},{\bm{Z}}_{0,c}\right]$ and $\left[{\bm{Z}}_{0}^{(l)},{\bm{Z}}_{0,c}^{(l)}\right]$ are orthonormal matrices. Below we will prove the induction hypotheses (28) in the base case when $t=0$ one by one.

The last term in (67) can be further bounded as

via Lemma 9. Plugging (68) into (67), we have

Setting $\delta=c\frac{\sigma_{r}^{3}\left({\bm{X}}^{\natural}\right)}{\sqrt{r}\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{3}}$ for a sufficiently small constant $c$, i.e. $m\gtrsim\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{6}}{\sigma_{r}^{6}\left({\bm{X}}^{\natural}\right)}nr\log n$, we get $\left\|{\bm{X}}_{0}{\bm{Q}}_{0}-{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}\allowbreak\leq C_{1}\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}$. Following similar procedures, we can also show $\left\|{\bm{X}}_{0}^{(l)}{\bm{Q}}_{0}^{(l)}-{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}\allowbreak\leq C_{1}\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}$.

G.2 Proof of (28b)

Following Weyl’s inequality, by (28a), we have

and similarly, $\left|\sigma_{i}\left({\bm{X}}_{0}^{(l)}\right)-\sigma_{i}\left({\bm{X}}^{\natural}\right)\right|\leq C_{1}\frac{\sigma_{r}^{2}\left({\bm{X}}^{\natural}\right)}{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}}$, for $i=1,\cdots,r$. Combined with Lemma 6, there exists some constant $c$ such that

We will bound each term in (69), respectively. For the first term, we have

where the last line follows from (66). Note that the first term in (70) can be bounded as

which follows Lemma 10 and Lemma 11. The second term in (70) can be bounded as

where (72) follows from Lemma 13 and the Davis-Kahan $\sin\Theta$ theorem , and (73) follows from Lemma 10 and Lemma 11.

where the first term of (74) is bounded similarly as (73), and (75) follows from Lemma 11. Combining (71), (73), and (75), we obtain

where the last inequality holds as long as $m\gtrsim\kappa\frac{\left\|{\bm{X}}^{\natural}\right\|_{{\mathsf{F}}}^{5}}{\sigma_{r}^{5}\left({\bm{X}}^{\natural}\right)}n\sqrt{r}\log{n}=O(nr^{3}\log n)$.

G.3 Proof of (28c)

with proper constants, following Lemma 37 in , we have

which implies that for every $1\leq l\leq m$ we can get

where (76) follows from Lemma 10 and Lemma 11, and (77) follows from (28b).

G.4 Finishing the Proof

The proof of Lemma 5 is now complete by appropriately adjusting the constants.

References