Fast Stochastic Algorithms for SVD and PCA: Convergence Properties and Convexity

Introduction

We consider the problem of recovering the top $k$ left singular vectors of a $d\times n$ matrix $X=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$, where $k\ll d$. This is equivalent to recovering the top $k$ eigenvectors of $XX^{\top}$, or equivalently, solving the optimization problem

This is one of the most fundamental matrix computation problems, and has numerous uses (such as low-rank matrix approximation and principal component analysis).

For large-scale matrices $X$, where exact eigendecomposition is infeasible, standard deterministic approaches are based on power iterations or variants thereof (e.g. the Lanczos method) . Alternatively, one can exploit the structure of Eq. (1) and apply stochastic iterative algorithms, where in each iteration we update a current $d\times k$ matrix $W$ based on one or more randomly-drawn columns $\mathbf{x}_{i}$ of $X$. Such algorithms have been known for several decades (), and enjoyed renewed interest in recent years, e.g. . Another stochastic approach is based on random projections, e.g. .

Unfortunately, each of these algorithms suffer from a different disadvantage: The deterministic algorithms are accurate (runtime logarithmic in the required accuracy $\epsilon$, under an eigengap condition), but require a full pass over the matrix for each iteration, and in the worst-case many such passes would be required (polynomial in the eigengap). On the other hand, each iteration of the stochastic algorithms is cheap, and their number is independent of the size of the matrix, but on the flip side, their noisy stochastic nature means they are not suitable for obtaining a high-accuracy solution (the runtime scales polynomially with $\epsilon$).

Recently, proposed a new practical algorithm, VR-PCA, for solving Eq. (1), which has a “best-of-both-worlds” property: The algorithm is based on cheap stochastic iterations, yet the algorithm’s runtime is logarithmic in the required accuracy $\epsilon$. More precisely, for the case $k=1$, $\mathbf{x}_{i}$ of bounded norm, and when there is an eigengap of $\lambda$ between the first and second leading eigenvalues of the covariance matrix $\frac{1}{n}XX^{\top}$, the required runtime was shown to be on the order of

The algorithm is therefore suitable for obtaining high accuracy solutions (the dependence on $\epsilon$ is logarithmic), but essentially at the cost of only $\mathcal{O}(\log(1/\epsilon))$ passes over the data. The algorithm is based on a recent variance-reduction technique designed to speed up stochastic algorithms for convex optimization problems (), although the optimization problem in Eq. (1) is inherently non-convex. See Section 3 for a more detailed description of this algorithm, and for more discussions as well as empirical results.

The results and analysis in left several issues open. For example, it is not clear if the quadratic dependence on $1/\lambda$ in Eq. (2) is necessary, since it is worse than the linear (or better) dependence that can be obtained with the deterministic algorithms mentioned earlier, as well as analogous results that can be obtained with similar techniques for convex optimization problems (where $\lambda$ is the strong convexity parameter). Also, the analysis was only shown for the case $k=1$, whereas often in practice, we may want to recover $k>1$ singular vectors simultaneously. Although proposed a variant of the algorithm for that case, and studied it empirically, no analysis was provided. Finally, the convergence guarantee assumed that the algorithm is initialized from a point closer to the optimum than what is attained with standard random initialization. Although one can use some other, existing stochastic algorithm to do this “warm-start”, no end-to-end analysis of the algorithm, starting from random initialization, was provided.

In this paper, we study these and related questions, and make the following contributions:

We propose a variant of VR-PCA to handle the $k>1$ case, and formally analyze its convergence (Section 3). The extension to $k>1$ is non-trivial, and requires tracking the evolution of the subspace spanned by the current solution at each iteration.

In Section 4, we study the convergence of VR-PCA starting from a random initialization. And show that with a slightly smarter initialization – essentially, random initialization followed by a single power iteration – the convergence results can be substantially improved. In fact, a similar initialization scheme should assist in the convergence of other stochastic algorithms for this problem, as long as a single power iteration can be performed.

In Section 5, we study whether functions similar to Eq. (1) have hidden convexity properties, which would allow applying existing convex optimization tools as-is, and improve the required runtime. For the $k=1$ case, we show that this is in fact true: Close enough to the optimum, and on a suitably-designed convex set, such a function is indeed $\lambda$-strongly convex. Unfortunately, the distance from the optimum has to be $\mathcal{O}(\lambda)$, and this precludes a better runtime in most practical regimes. However, it still indicates that a better runtime and dependence on $\lambda$ should be possible.

Some Preliminaries and Notation

We consider a $d\times n$ matrix $X$ composed of $n$ columns $(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$, and let

Thus, Eq. (1) is equivalent to finding the $k$ leading eigenvectors of $A$.

The VR-PCA Algorithm and a Block Version

We begin by recalling the algorithm of for the $k=1$ case (Algorithm 1), and then discuss its generalization for $k>1$.

To handle the $k>1$ case (where more than one eigenvector should be recovered), one simple technique is deflation, where we recover the leading eigenvectors $\mathbf{v}_{1},\mathbf{v}_{2},\ldots,\mathbf{v}_{k}$ one-by-one, each time using the $k=1$ algorithm. However, a disadvantage of this approach is that it requires a positive eigengap between all top $k$ eigenvalues, otherwise the algorithm is not guaranteed to converge. Thus, an algorithm which simultaneously recovers all $k$ leading eigenvectors is preferable.

We now turn to provide a formal analysis of Algorithm 2, which directly generalizes the analysis of Algorithm 1 given in :

Define the $d\times d$ matrix $A$ as $\frac{1}{n}XX^{\top}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}$, and let $V_{k}$ denote the $d\times k$ matrix composed of the eigenvectors corresponding to the largest $k$ eigenvalues. Suppose that

$\max_{i}\|\mathbf{x}_{i}\|^{2}\leq r$ for some $r>0$.

$A$ has eigenvalues $s_{1}>s_{2}\geq\ldots\geq s_{d}$, where $s_{k}-s_{k+1}=\lambda$ for some $\lambda>0$.

Let $\delta,\epsilon\in(0,1)$ be fixed. If we run the algorithm with any epoch length parameter $m$ and step size $\eta$, such that

(where $c,c^{\prime},c^{\prime\prime}$ designate certain positive numerical constants), and for $T=\left\lceil\frac{\log(1/\epsilon)}{\log(2/\delta)}\right\rceil$ epochs, then with probability at least $1-\lceil\log_{2}(1/\epsilon)\rceil\delta$, it holds that

For any orthogonal $W$, $k-\|V_{k}^{\top}W\|_{F}^{2}$ lies between and $k$, and equals when the column spaces of $V_{k}$ and $W$ are the same (i.e., when $W$ spans the $k$ leading singular vectors). According to the theorem, taking appropriateSpecifically, we can take $m=c^{\prime}\log(2/\delta)/\eta\lambda$ and $\eta=a\delta^{2}/r^{2}\lambda$, where $a$ is sufficiently small to ensure that the first and third condition in Eq. (3) holds. It can be verified that it’s enough to take $a=\min\left\{c,\frac{c^{\prime\prime}}{4\delta^{2}ck\log(2/\delta)},\frac{1}{4\delta^{2}c}\left(\frac{c^{\prime\prime}}{k\log(2/\delta)}\right)^{2}\right\}$. $\eta=\Theta(\lambda/(kr)^{2})$, and $m=\Theta((rk/\lambda)^{2})$, the algorithm converges with high probability to a high-accuracy approximation of $V_{k}$. Moreover, the runtime of each epoch of the algorithm equals $\mathcal{O}(mdk^{2}+dnk)$. Overall, we get the following corollary:

This runtime bound is the same showed that it’s possible to further improve the runtime for sparse $X$, replacing $d$ by the average column sparsity $d_{s}$. This is done by maintaining parameters in an implicit form, but it’s not clear how to implement a similar trick in the block version, where $k>1$. as that of for $k=1$.

Warm-Start and the Power of a Power Iteration

In this section, we study the runtime required to compute a starting point satisfying the conditions of Theorem 1, starting from a random initialization. Combined with Theorem 1, this gives us an end-to-end analysis of the runtime required to find an $\epsilon$-accurate solution, starting from a random point. For simplicity, we will only discuss the case $k=1$, i.e. where our goal is to compute the single leading eigenvector $\mathbf{v}_{1}$, although our observations can be generalized to $k>1$. In the $k=1$ case, Theorem 1 kicks in once we find a vector $\mathbf{w}$ satisfying $\langle\mathbf{v}_{1},\mathbf{w}\rangle^{2}\geq\frac{1}{2}$.

To the best of our knowledge, the existing iteration complexity guarantees for such algorithms (assuming the norm constraint $r\leq 1$ for simplicity) scale at leastFor example, this holds for , although the bound only guarantees the existence of some iteration which produces the desired output. The guarantee of scale as $d^{2}/\lambda^{2}$, and the guarantee of scales as $d/\lambda^{3}$ in our setting. as $d/\lambda^{2}$. Since the runtime of each iteration is $\mathcal{O}(d)$, we get an overall runtime of $\mathcal{O}((d/\lambda)^{2})$.

The dependence on $d$ in the iteration bound stems from the fact that with a random initial unit vector $\mathbf{w}_{0}$, we have $\langle\mathbf{v}_{1},\mathbf{w}_{0}\rangle^{2}\approx\frac{1}{d}$. Thus, we begin with a vector almost orthogonal to the leading eigenvector $\mathbf{v}_{1}$ (depending on $d$). In a purely stochastic setting, where only noisy information is available, this necessitates conservative updates at first, and in all the analyses we are aware of, the number of iterations appear to necessarily scale at least linearly with $d$.

For $\mathbf{w}_{0}$ as above, it holds for any $\delta$ that with probability at least $1-\frac{1}{d}-\delta$,

where $\text{nrank}(A)=\frac{\|A\|_{F}^{2}}{\|A\|_{sp}^{2}}$ is the numerical rank of $A$.

The numerical rank (see e.g. ) is a relaxation of the standard notion of rank: For any $d\times d$ matrix $A$, nrank(A) is at most the rank of $A$ (which in turn is at most $d$). However, it will be small even if $A$ is just close to being low-rank. In many if not most machine learning applications, we are interested in matrices which tend to be approximately low-rank, in which case nrank(A) is much smaller than $d$ or even a constant. Therefore, by a single power iteration, we get an initial point $\mathbf{w}_{0}$ for which $\langle\mathbf{v}_{1},\mathbf{w}_{0}\rangle^{2}$ is on the order of $1/\text{nrank(A)}$, which can be much larger than the $1/d$ given by a random initialization, and is never substantially worse.

Let $s_{1}\geq s_{2}\geq\ldots\geq s_{d}\geq 0$ be the $d$ eigenvalues of $A$, with eigenvectors $\mathbf{v}_{1},\ldots,\mathbf{v}_{d}$. We have

Since $\mathbf{w}$ is distributed according to a standard Gaussian distribution, which is rotationally symmetric, we can assume without loss of generality that $\mathbf{v}_{1},\ldots,\mathbf{v}_{d}$ correspond to the standard basis vectors $\mathbf{e}_{1},\ldots,\mathbf{e}_{d}$, in which case the above reduces to

where $w_{1},\ldots,w_{d}$ are independent and scalar random variables with a standard Gaussian distribution.

First, we note that $s_{1}^{2}$ equals $\|A\|_{sp}^{2}$, the spectral norm of $A$, whereas $\sum_{i=1}^{d}s_{i}^{2}$ equals $\|A\|_{F}^{2}$, the Frobenius norm of $A$. Therefore, $\frac{s_{1}^{2}}{\sum_{i}s_{i}^{2}}=\frac{\|A\|_{sp}^{2}}{\|A\|_{F}^{2}}=\frac{1}{\text{nrank}(A)}$, and we get overall that

We consider the random quantity $w_{1}^{2}/\max_{i}w_{i}^{2}$, and independently bound the deviation probability of the numerator and denominator. First, for any $t\geq 0$ we have

Combining Eq. (5) and Eq. (6), with a union bound, we get that for any $t_{1},t_{2}\geq 0$, it holds with probability at least $1-\sqrt{\frac{2}{\pi}t_{1}}-\exp(-t_{2}^{2}/2)$ that

To slightly simplify this for readability, we take $t_{2}=\sqrt{2\log(d)}$, and substitute $\delta=\sqrt{\frac{2}{\pi}t_{1}}$. This implies that with probability at least $1-\delta-1/d$,

Plugging back into Eq. (4), the result follows. ∎

This result can be plugged into the existing analyses of purely stochastic PCA/SVD algorithms, and can often improve the dependence on the $d$ factor in the iteration complexity bounds to a dependence on the numerical rank of $A$. We again emphasize that this is applicable in a situation where we can actually perform a power iteration, and not in a purely stochastic setting where we only have access to an i.i.d. data stream (nevertheless, it would be interesting to explore whether this idea can be utilized in such a streaming setting as well).

To give a concrete example of this, we provide a convergence analysis of the VR-PCA algorithm (Algorithm 1), starting from an arbitrary initial point, bounding the total number of stochastic iterations required by the algorithm in order to produce a point satisfying the conditions of Theorem 1 (from which point the analysis of Theorem 1 takes over). Combined with Theorem 1, this analysis also justifies that VR-PCA indeed converges starting from a random initialization.

(for some universal constant $c$). Then with probability at least $1-\delta$, after

stochastic iterations (lines $6-10$ in the pseudocode, where $c^{\prime}$ is again a universal constant), we get a point $\mathbf{w}_{T}$ satisfying $1-\langle\mathbf{v}_{1},\mathbf{w}_{T}\rangle^{2}\leq\frac{1}{2}$. Moreover, if $\eta$ is chosen on the same order as the upper bound in Eq. (7), then

Note that the analysis does not depend on the choice of the epoch size $m$, and does not use the special structure of VR-PCA (in fact, the technique we use is applicable to any algorithm which takes stochastic gradient steps to solve this type of problemAlthough there exist previous analyses of such algorithms in the literature, they unfortunately do not quite apply to our algorithm, for various technical reasons.). The proof of the theorem appears in Section 6.2.

By Corollary 1, the runtime required by VR-PCA from that point to get an $\epsilon$-accurate solution is

so the sum of the two expressions (which is $d\left(n+\frac{1}{\lambda^{2}}\right)$ up to log-factors), represents the total runtime required by the algorithm.

Convexity and Non-Convexity of the Rayleigh Quotient

As mentioned in the introduction, an intriguing open question is whether the $d\left(n+\frac{1}{\lambda^{2}}\right)\log\left(\frac{1}{\epsilon}\right)$ runtime guarantees from the previous sections can be further improved. Although a linear dependence on $d,n$ seems unavoidable, this is not the case for the quadratic dependence on $1/\lambda$. Indeed, when using deterministic methods such as power iterations or the Lanczos method, the dependence on $\lambda$ in the runtime is only $1/\lambda$ or even $\sqrt{1/\lambda}$ . In the world of convex optimization from which our algorithmic techniques are derived, the analog of $\lambda$ is the strong convexity parameter of the function, and again, it is possible to get a dependence of $1/\lambda$, or even $\sqrt{1/\lambda}$ with accelerated schemes (see e.g. in the context of the variance-reduction technique we use). Is it possible to get such a dependence for our problem as well?

Another question is whether the non-convex problem that we are tackling (Eq. (1)) is really that non-convex. Clearly, it has a nice structure (since we can solve the problem in polynomial time), but perhaps it actually has hidden convexity properties, at least close enough to the optimal points? We note that Eq. (1) can be “trivially” convexified, by re-casting it as an equivalent semidefinite program . However, that would require optimization over $d\times d$ matrices, leading to poor runtime and memory requirements. The question here is whether we have any convexity with respect to the original optimization problem over “thin” $d\times k$ matrices.

In fact, the two questions of improved runtime and convexity are closely related: If we can show that the optimization problem is convex in some domain containing an optimal point, then we may be able to use fast stochastic algorithms designed for convex optimization problems, inheriting their good guarantees.

To discuss these questions, we will focus on the $k=1$ case for simplicity (i.e., our goal is to find a leading eigenvector of the matrix $A=\frac{1}{n}XX^{\top}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i}\mathbf{x}_{i}^{\top}$), and study potential convexity properties of the negative Rayleigh quotient,

Note that for $k=1$, this function coincides with Eq. (1) on the unit Euclidean sphere, and with the same optimal points, but has the nice property of being defined on the entire Euclidean space (thus, at least its domain is convex).

At a first glance, such functions $F_{A}$ appear to potentially be convex at some bounded distance from an optimum, as illustrated for instance in the case where A=\left(\begin{array}[]{cc}1&0\\ 0&0\\ \end{array}\right) (see Figure 1). Unfortunately, it turns out that the figure is misleading, and in fact the function is not convex almost everywhere:

For the matrix $A$ above, the Hessian of $F_{A}$ is not positive semidefinite for all but a measure-zero set.

The leading eigenvector of $A$ is $\mathbf{v}_{1}=(1,0)$, and $F_{A}(\mathbf{w})=-\frac{w_{1}^{2}}{w_{1}^{2}+w_{2}^{2}}$. The Hessian of this function at some $\mathbf{w}$ equals

The determinant of this $2\times 2$ matrix equals

The theorem implies that we indeed cannot use convex optimization tools as-is on the function $F_{A}$, even if we’re close to an optimum. However, the non-convexity was shown for $F_{A}$ as a function over the entire Euclidean space, so the result does not preclude the possibility of having convexity on a more constrained, lower-dimensional set. In fact, this is what we are going to do next: We will show that if we are given some point $\mathbf{w}_{0}$ close enough to an optimum, then we can explicitly construct a simple convex set, such that

The set includes an optimal point of $F_{A}$.

The function $F_{A}$ is $\mathcal{O}(1)$-smooth and $\lambda$-strongly convex in that set.

This means that we can potentially use a two-stage approach: First, we use some existing algorithm (such as VR-PCA) to find $\mathbf{w}_{0}$, and then switch to a convex optimization algorithm designed to handle functions with a finite sum structure (such as $F_{A}$). Since the runtime of such algorithms scale better than VR-PCA, in terms of the dependence on $\lambda$, we can hope for an overall runtime improvement.

Unfortunately, this has a catch: To make it work, we need to have $\mathbf{w}_{0}$ very close to the optimum – in fact, we require $\|\mathbf{v}_{1}-\mathbf{w}_{0}\|\leq\mathcal{O}(\lambda)$, and we show (in Theorem 5) that such a dependence on the eigengap $\lambda$ cannot be avoided (perhaps up to a small polynomial factor). The issue is that the runtime to get such a $\mathbf{w}_{0}$, using stochastic-based approaches we are aware of, would scale at least quadratically with $1/\lambda$, but getting dependence better than quadratic was our problem to begin with. For example, the runtime guarantee using VR-PCA to get such a point $\mathbf{w}_{0}$ (even if we start from a good point as specified in Theorem 1) is on the order of

whereas the best known guarantees on getting an $\epsilon$-optimal solution for $\lambda$-strongly convex and smooth functions (see ) is on the order of

Therefore, the total runtime we can hope for would be on the order of

In comparison, the runtime guarantee of using just VR-PCA to get an $\epsilon$-accurate solution is on the order of

Unfortunately, Eq. (9) is the same as Eq. (8) up to log-factors, and the difference is not significant unless the required accuracy $\epsilon$ is extremely small (exponentially small in $n,1/\lambda$). Therefore, our construction is mostly of theoretical interest. However, it still shows that asymptotically, as $\epsilon\rightarrow 0$, it is indeed possible to have runtime scaling better than Eq. (9). This might hint that designing practical algorithms, with better runtime guarantees for our problem, may indeed be possible.

To explain our construction, we need to consider two convex sets: Given a unit vector $\mathbf{w}_{0}$, define the hyperplane tangent to $\mathbf{w}_{0}$,

as well as a Euclidean ball of radius $r$ centered at $\mathbf{w}_{0}$:

The convex set we use, given such a $\mathbf{w}_{0}$, is simply the intersection of the two, $H_{\mathbf{w}_{0}}\cap B_{\mathbf{w}_{0}}(r)$, where $r$ is a sufficiently small number (see Figure 2).

The following theorem shows that if $\mathbf{w}_{0}$ is $\mathcal{O}(\lambda)$-close to an optimal point (a leading eigenvector $\mathbf{v}_{1}$ of $A$), and we choose the radius of $B_{w_{0}}(r)$ appropriately, then $H_{\mathbf{w}_{0}}\cap B_{\mathbf{w}_{0}}(r)$ contains an optimal point, and the function $F_{A}$ is indeed $\lambda$-strongly convex and smooth on that set. For simplicity, we will assume that $A$ is scaled to have spectral norm of $1$, but the result can be easily generalized.

For any positive semidefinite $A$ with spectral norm $1$, eigengap $\lambda$ and a leading eigenvector $\mathbf{v}_{1}$, and any unit vector $\mathbf{w}_{0}$ such that $\|\mathbf{w}_{0}-\mathbf{v}_{1}\|\leq\frac{\lambda}{44}$, the function $F_{A}(\mathbf{w})$ is $20$-smooth and $\lambda$-strongly convex on the convex set $H_{\mathbf{w}_{0}}\cap B_{\mathbf{w}_{0}}\left(\frac{\lambda}{22}\right)$, which contains a global optimum of $F_{A}$.

The proof of the theorem appears in Subsection 6.3. Finally, we show below that a polynomial dependence on the eigengap $\lambda$ is unavoidable, in the sense that the convexity property is lost if $\mathbf{w}_{0}$ is significantly further away from $\mathbf{v}_{1}$.

For any $\lambda,\epsilon\in\left(0,\frac{1}{2}\right)$, there exists a positive semidefinite matrix $A$ with spectral norm $1$, eigengap $\lambda$, and leading eigenvector $\mathbf{v}_{1}$, as well as a unit vector $\mathbf{w}_{0}$ for which $\|\mathbf{v}_{1}-\mathbf{w}_{0}\|\leq\sqrt{2(1+\epsilon)\lambda)}$, such that $F_{A}$ is not convex in any neighborhood of $\mathbf{w}_{0}$ on $H_{\mathbf{w}_{0}}$.

for which $\mathbf{v}_{1}=(1,0,0)$, and take

where $p=\sqrt{(1+\epsilon)\lambda}$ (which ensures $\|\mathbf{v}_{1}-\mathbf{w}_{0}\|^{2}=\sqrt{2p^{2}}=\sqrt{2(1+\epsilon)\lambda}$). Consider the ray $\{(\sqrt{1-p^{2}},t,p):t\geq 0\}$, and note that it starts from $\mathbf{w}_{0}$ and lies in $H_{\mathbf{w}_{0}}$. The function $F_{A}$ along that ray (considering it as a function of $t$) is of the form

The second derivative with respect to $t$ equals

where we plugged in the definition of $p$. This is a negative quantity for any $t<\frac{1}{\sqrt{3}}$. Therefore, the function $F_{A}$ is strictly concave (and not convex) along the ray we have defined and close enough to $\mathbf{w}_{0}$, and therefore isn’t convex in any neighborhood of $\mathbf{w}_{0}$ on $H_{\mathbf{w}_{0}}$. ∎

Proofs

Although the proof structure generally mimics the proof of Theorem 1 in for the $k=1$ special case, it is more intricate and requires several new technical tools. To streamline the presentation of the proof, we begin with proving a series of auxiliary lemmas in Subsection 6.1.1, and then move to the main proof in Subsection 6.1. The main proof itself is divided into several steps, each constituting one or more lemmas.

For any $B,C,D\succeq 0$, it holds that $\operatorname{Tr}(BC)\geq\operatorname{Tr}(B(C-D))$ and $\operatorname{Tr}(BC)\geq\operatorname{Tr}((B-D)C)$.

It is enough to prove that for any positive semidefinite matrices $E,G$, it holds that $\operatorname{Tr}(EG)\geq 0$. The lemma follows by taking either $E=B,G=D$ (in which case, $\operatorname{Tr}(BC)=\operatorname{Tr}(B(C-D))+\operatorname{Tr}(BD)\geq\operatorname{Tr}(B(C-D))$), or $E=D,G=C$ (in which case, $\operatorname{Tr}(BC)=\operatorname{Tr}((B-D)C)+\operatorname{Tr}(DC)\geq\operatorname{Tr}((B-D)C)$).

Any positive semidefinite matrix $M$ can be written as the product $M^{1/2}M^{1/2}$ for some symmetric matrix $M^{1/2}$ (known as the matrix square root of $M$). Therefore,

We begin by proving the one-dimensional case, where $B,C$ are scalars $b\geq 0,c>0$. The inequality then becomes $bc^{-1}\geq b(2-c)$, which is equivalent to $1\geq c(2-c)$, or upon rearranging, $(c-1)^{2}\geq 0$, which trivially holds.

Turning to the general case, we note that by Lemma 2, it is enough to prove that $C^{-1}-(2I-C)\succeq 0$. To prove this, we make a couple of observations. The positive definite matrix $C$ (like any positive definite matrix) has a singular value decomposition which can be written as $USU^{\top}$, where $U$ is an orthogonal matrix, and $S$ is a diagonal matrix with positive entries. Its inverse is $US^{-1}U^{\top}$, and $2I-C=2I-USU^{\top}=U(2I-S)U^{\top}$. Therefore,

To show this matrix is positive semidefinite, it is enough to show that each diagonal entry of $S^{-1}-(2I-S)$ is non-negative. But this reduces to the one-dimensional result we already proved, when $b=1$ and $c>0$ is any diagonal entry in $S$. Therefore, $C^{-1}-(2I-C)\succeq 0$, from which the result follows. ∎

The first inequality is immediate from Cauchy-Shwartz. As to the second inequality, letting $\mathbf{c}_{i}$ denote the $i$-th column of $C$, and $\|\cdot\|_{2}$ the Euclidean norm for vectors,

Let $B_{1},B_{2},Z_{1},Z_{2}$ be $k\times k$ square matrices, where $B_{1},B_{2}$ are fixed and $Z_{1},Z_{2}$ are stochastic and zero-mean (i.e. their expectation is the all-zeros matrix). Furthermore, suppose that for some fixed $\alpha,\gamma,\delta>0$, it holds with probability $1$ that

For all $\nu\in$, $B_{2}+\nu Z_{2}\succeq\delta I$.

$\max\{\|Z_{1}\|_{F},\|Z_{2}\|_{F}\}\leq\alpha$.

Since $B_{2}+\nu Z_{2}$ is positive definite, it is always invertible, hence $f(\nu)$ is indeed well-defined. Moreover, it can be differentiated with respect to $\nu$, and we have

Again differentiating with respect to $\nu$, we have

Using Lemma 4 and the triangle inequality, this is at most

Applying a Taylor expansion to $f(\cdot)$ around $\nu=0$, with a Lagrangian remainder term, and substituting the values for $f^{\prime}(\nu),f^{\prime\prime}(\nu)$, we can lower bound $f(1)$ as follows:

Taking expectation over $Z_{1},Z_{2}$, and recalling they are zero-mean, we get that

Let $U_{1},\ldots,U_{k}$ and $R_{1},R_{2}$ be positive semidefinite matrices, such that $R_{2}-R_{1}\succeq 0$, and define the function

over all $(x_{1}\ldots x_{k})\in[\alpha,\beta]^{d}$ for some $\beta\geq\alpha\geq 0$. Then $~{}\min_{(x_{1}\ldots x_{k})\in[\alpha,\beta]^{d}}f(\mathbf{x})=f(\alpha,\ldots,\alpha)$.

Taking a partial derivative of $f$ with respect to some $x_{j}$, we have

By the lemma’s assumptions, each matrix in the product above is positive semidefinite, hence the product is positive semidefinite, and the trace is non-negative. Therefore, $\frac{\partial}{\partial x_{j}}f(\mathbf{x})\geq 0$, which implies that the function is minimized when each $x_{j}$ takes its smallest possible value, i.e. $\alpha$. ∎

Let $B$ be a $k\times k$ matrix with minimal singular value $\delta$. Then

so it remains to prove $1-\frac{\|B^{\top}B\|_{F}^{2}}{\|B\|_{F}^{2}}\geq\frac{\delta^{2}}{k}\left(k-\|B\|_{F}^{2}\right)$. Let $\sigma_{1},\ldots,\sigma_{k}$ denote the vector of singular values of $B$. The singular values of $B^{\top}B$ are $\sigma_{1}^{2},\ldots,\sigma_{k}^{2}$, and the Frobenius norm of a matrix equals the Euclidean norm of its vector of singular values. Therefore, the lemma is equivalent to requiring

assuming $\sigma_{i}\in[\delta,1]$ for all $i$. This holds since

For any $d\times k$ matrices $C,D$ with orthonormal columns, let

be the nearest orthonormal-columns matrix to $C$ in the column space of $D$ (where $B$ is a $k\times k$ matrix). Then the matrix $B$ minimizing the above equals $B=VU^{\top}$, where $C^{\top}D=USV^{\top}$ is the SVD decomposition of $C^{\top}D$, and it holds that

Since $D$ has orthonormal columns, we have $D^{\top}D=I$, so the definition of $B$ is equivalent to

This is the orthogonal Procrustes problem (see e.g. ), and the solution is easily shown to be $B=VU^{\top}$ where $USV^{\top}$ is the SVD decomposition of $C^{\top}D$. In this case, and using the fact that $\|C\|_{F}^{2}=\|D\|_{F}^{2}=k$ (as $C,D$ have orthonormal columns), we have that $\|C-D_{C}\|_{F}^{2}$ equals

Since the trace function is similarity-invariant, this equals $2k-\operatorname{Tr}(S)$. Let $s_{1}\ldots,s_{k}$ be the diagonal elements of $S$, and note that they can be at most $1$ (since they are the singular values of $C^{\top}D$, and both $C$ and $D$ have orthonormal columns). Recalling that the Frobenius norm equals the Euclidean norm of the singular values, we can therefore upper bound the above as follows:

Let $W_{t},W^{\prime}_{t}$ be as defined in Algorithm 2, where we assume $\eta<\frac{1}{3}$. Then for any $d\times k$ matrix $V_{k}$ with orthonormal columns, it holds that

Letting $\mathbf{s}_{t},\mathbf{s}_{t-1}$ denote the vectors of singular values of $V_{k}^{\top}W_{t}$ and $V_{k}^{\top}W_{t-1}$, and noting that they are both in $^{k}$ (as $V_{k},W_{t-1},W_{t}$ all have orthonormal columns), the left hand side of the inequality in the lemma statement equals

where $\|\cdot\|_{\infty}$ is the infinity norm. By Weyl’s matrix perturbation theoremUsing its version for singular values, which implies that the singular values of matrices $B$ and $B+E$ are different by at most $\|E\|_{sp}$. , this is upper bounded by

Recalling the relationship between $W_{t}$ and $W_{t-1}$ from Algorithm 2, we have that

Plugging back to Eq. (10), the result follows. ∎

1.2 Main Proof

To simplify the technical derivations, note that the algorithm remains the same if we divide each $\mathbf{x}_{i}$ by $\sqrt{r}$, and multiply $\eta$ by $r$. Since $\max_{i}\|\mathbf{x}_{i}\|^{2}\leq r$, this corresponds to running the algorithm with step-size $\eta r$ rather than $\eta$, on a re-scaled dataset of points with squared norm at most $1$, and with an eigengap of $\lambda/r$ instead of $\lambda$. Therefore, we can simply analyze the algorithm assuming that $\max_{i}\|\mathbf{x}_{i}\|^{2}\leq 1$, and in the end plug in $\lambda/r$ instead of $\lambda$, and $\eta r$ instead of $\eta$, to get a result which holds for data with squared norm at most $r$.

Part I: Establishing a Stochastic Recurrence Relation

We begin by focusing on a single iteration $t$ of the algorithm, and analyze how $\|V_{k}^{\top}W_{t}\|_{F}^{2}$ (which measures the similarity between the column spaces of $V_{k}$ and $W_{t}$) evolves during that iteration. The key result we need is Lemma 10 below, which is specialized for our algorithm in Lemma 11.

Let $A$ be a $d\times d$ symmetric matrix with all eigenvalues $s_{1}\geq s_{2}\geq\ldots\geq s_{d}$ in $$, and suppose that$s_{k}-s_{k+1}\geq\lambda$for some$\lambda>0$.

Let $N$ be a $d\times k$ zero-mean random matrix such that $\|N\|_{F}\leq\sigma^{F}_{N}$ and $\|N\|_{sp}\leq\sigma^{sp}_{N}$ with probability $1$, and define

Let $W$ be a $d\times k$ matrix with orthonormal columns, and define

for some $\eta\in\left[0,\frac{1}{4\max\{1,\sigma^{F}_{N}\}}\right]$.

If $V_{k}=[\mathbf{v}_{1},\mathbf{v}_{2}\ldots,\mathbf{v}_{k}]$ is the $d\times k$ matrix of $A$’s first $k$ eigenvectors, then the following holds:

If $\|V_{k}^{\top}W\|_{F}^{2}\geq k-\frac{1}{2}$, then

Using the fact that $\operatorname{Tr}(BCD)=\operatorname{Tr}(CDB)$ for any matrices $B,C,D$, we have

$Z_{1},Z_{2}$ are zero mean: This holds since they are linear in $N$, and $N$ is assumed to be zero-mean.

$B_{2}+\nu Z_{2}\succeq\frac{3}{8}I$ for all $\nu\in$: Recalling the definition of $B_{2},Z_{2}$, and the facts that $A\succeq 0$, $N^{\top}N\succeq 0$ (by construction), and $W^{\top}W=I$, we have that $B_{2}\succeq I$. Moreover, the spectral norm of $Z_{2}$ is at most

which by the assumption on $\eta$ is at most $2\frac{1}{4}\left(1+\frac{1}{4}\right)=\frac{5}{8}$. This implies that the smallest singular value of $B_{2}+\nu Z_{2}$ is at least $1-\nu(5/8)\geq 3/8$.

$\max\{\|Z_{1}\|_{F},\|Z_{2}\|_{F}\}\leq\frac{5}{2}\eta\sigma^{F}_{N}$: By definition of $Z_{1},Z_{2}$, and using Lemma 4, the Frobenius norm of these two matrices is at most

which by the assumption on $\eta$ is at most $2\eta\sigma^{F}_{N}\left(1+\frac{1}{4}\right)=\frac{5}{2}\eta\sigma^{F}_{N}$.

$\|B_{1}+\eta Z_{1}\|_{sp}\leq\left(\frac{1}{4}\sigma^{sp}_{N}+2\right)^{2}$: Using the definition of $B_{1},Z_{1}$ and the assumption $\eta\leq\frac{1}{4}$,

Applying Lemma 5 and plugging back to Eq. (11), we get

We now turn to lower bound $\operatorname{Tr}\left(B_{1}B_{2}^{-1}\right)$, by first re-writing $B_{1},B_{2}$ in a different form. For $i=1,\ldots,d$, let

where $\mathbf{v}_{i}$ is the eigenvector of $A$ corresponding to the eigenvalue $s_{i}$. Note that each $U_{i}$ is positive semidefinite, and $\sum_{i=1}^{d}U_{i}=W^{\top}W=I$. We have

Plugging Eq. (13) and Eq. (14) back into Eq. (12), we get

Recalling that $s_{1}\geq s_{2}\geq\ldots\geq s_{k}$ and letting $\alpha=(1+\eta s_{k})^{2},\beta=(1+\eta s_{1})^{2}$, the trace term can be lower bounded by

Applying Lemma 6 (noting that as required by the lemma, $\sum_{i=k+1}^{d}(1+\eta s_{i})^{2}U_{i}+\eta^{2}N^{\top}N-\eta^{2}N^{\top}V_{k}V_{k}^{\top}N=\sum_{i=k+1}^{d}(1+\eta s_{i})^{2}U_{i}+\eta^{2}N^{\top}\left(I-V_{k}V_{k}^{\top}\right)N\succeq 0$), we can lower bound the above by

Using Lemma 2, this can be lower bounded by

Recalling that $I=\sum_{i=1}^{d}U_{i}=\sum_{i=1}^{k}U_{i}+\sum_{i=k+1}^{d}U_{i}$, this can be simplified to

Since $U_{i}\succeq 0$, then using Lemma 3, we can lower bound the expression above by shrinking each of the $\left(2-\left(\frac{1+\eta s_{i}}{1+\eta s_{k}}\right)^{2}\right)$ terms. In particular, since $s_{i}\leq s_{k}-\lambda$ for each $i\geq k+1$,

which by the assumption that $\eta\leq 1/4$ and $s_{k}\leq s_{1}\leq 1$, is at least $1+\frac{4}{5}\eta\lambda$. Plugging this back into Eq. (16), and recalling that $\sum_{i=1}^{d}U_{i}=I$, we get the lower bound

Recall that this is a lower bound on the trace term in Eq. (15). Plugging it back and slightly simplifying, we get

The trace term above can be re-written (using the definition of $U_{i}$ and the fact that $\operatorname{Tr}(B^{\top}B)=\|B\|_{F}^{2}$) as

Applying Lemma 7, and letting $\delta$ denote the minimal singular value of $V_{k}^{\top}W$, this is lower bounded by

Taking the first argument of the max term in Eq. (17), we get

Subtracting $1$ from both sides and simplifying, we get

Suppose that $\|V_{k}^{\top}W\|_{F}^{2}\geq k-\frac{1}{2}$. Taking the second argument of the max term in Eq. (17), we get

Subtracting both sides from $k$, , we get

Since $k\geq 1$, we can lower bound the $\left(k-\frac{1}{2}\right)$ term by $\frac{k}{2}$. Moreover, the condition $k-\|V_{k}^{\top}W\|_{F}^{2}\leq\frac{1}{2}$ implies that the singular values $\sigma_{1},\ldots,\sigma_{k}$ of $V_{k}^{\top}W$ satisfy $k-\sum_{i=1}^{k}\sigma_{i}^{2}\leq\frac{1}{2}$. But each $\sigma_{i}$ is in $$(as$V_{k},W$have orthonormal columns), so no$\sigma_{i}$can be less than$\frac{1}{2}$. This implies that$\delta\geq\frac{1}{2}$. Plugging the lower bounds$k-\frac{1}{2}\geq\frac{k}{2}$and$\delta\geq\frac{1}{2}$ into the above, we get

Let $A,W_{t}$ be as defined in Algorithm 2, and suppose that $\eta\in\left[0,\frac{1}{23\sqrt{k}}\right]$. Then the following holds for some positive numerical constants $c_{1},c_{2},c_{3}$:

If $\|V_{k}^{\top}W_{t}\|_{F}^{2}\geq k-\frac{1}{2}$, then

so we may take $\sigma^{sp}_{N}=4$. As to the Frobenius norm, using Lemma 4 and a similar calculation, we have

to be the nearest orthonormal-columns matrix to $W_{t}$ in the column space of $V_{k}$, and

Plugging $\sigma^{sp}_{N}$ and $(\sigma^{F}_{N})^{2}$ into the $r_{N}$ as defined in Lemma 10, and picking any $\eta\in[0,\frac{1}{23\sqrt{k}}]$ (which satisfies the condition in Lemma 10 that $\eta\in\left[0,\frac{1}{4\max\{1,\sigma^{F}_{N}\}}\right]$, since $4\max\{1,\sigma^{F}_{n}\}\leq 4\max\{1,\sqrt{16*2k}\}<23\sqrt{k}$), we get

This implies that $r_{N}\leq 36800k$ always, which by application of Lemma 10, gives the first part of our lemma. As to the second part, assuming $\|V_{k}^{\top}W_{t}\|_{F}^{2}\geq k-\frac{1}{2}$ and applying Lemma 10, we get that

This corresponds to the lemma statement. ∎

Part II: Solving the Recurrence Relation for a Single Epoch

Suppose that $\eta=\alpha\lambda$, where $\alpha$ is a sufficiently small constant to be chosen later. Also, let

Then Lemma 11 tells us that if $\alpha$ is a sufficiently small constant, $b_{t}\leq\frac{1}{2}$, then

for some numerical constants $c,c^{\prime}$.

Let $B$ be the event that $b_{t}\leq\frac{1}{2}$ for all $t=0,1,2,\ldots,m$. Then for certain positive numerical constants $c_{1},c_{2},c_{3}$, if $\alpha\leq c_{1}$, then

where the expectation is over the randomness in the current epoch.

Note that the first equality holds, since conditioned on $W_{t}$, $b_{t+1}$ is independent of $b_{1},\ldots,b_{t}$, so the event $B$ is equivalent to just requiring $b_{t+1}\leq 1/2$.

Taking expectation over $W_{t}$ (conditioned on $B$), we get that

We now turn to prove that the event $B$ assumed in Lemma 12 indeed holds with high probability:

The following holds for certain positive numerical constants $c_{1},c_{2},c_{3}$: If $\alpha\leq c_{1}$, then for any $\beta\in(0,1)$ and $m$, if

then it holds with probability at least $1-\beta$ that

$|b_{t+1}-b_{t}|$ is bounded by $c^{\prime}_{3}k\alpha\lambda$ for some constant $c^{\prime}_{3}$: Applying Lemma 9, and assuming that $\alpha$ is at most some sufficiently small constant $c_{1}$ (e.g. $\alpha\leq\frac{1}{12}$, so $\eta=\alpha\lambda\leq\frac{1}{12}$),

Armed with these facts, and using the maximal version of the Hoeffding-Azuma inequality , it follows that with probability at least $1-\beta$, it holds simultaneously for all $t=1,\ldots,m$ (and for $t=0$ by assumption) that

Combining Lemma 12 and Lemma 13, and using Markov’s inequality, we get the following corollary:

Let confidence parameters $\beta,\gamma\in(0,1)$ be fixed. Suppose that $m,\alpha$ are chosen such that $\alpha\leq c_{1}$ and

where $c_{1},c_{2},c_{3}$ are certain positive numerical constants. Then with probability at least $1-(\beta+\gamma)$, it holds that

for some positive numerical constants $c,c^{\prime}$.

Part III: Analyzing the Entire Algorithm’s Run

Given the analysis in Lemma 14 for a single epoch, we are now ready to prove our theorem. Let

then we get with probability at least $1-(\beta+\gamma)$ that

Using the inequality $(1-(1/x))^{ax}\leq\exp(-a)$, which holds for any $x>1$ and any $a$, and taking $x=1/(c\alpha\lambda^{2})$ and $a=3\log(1/\gamma)$, we can upper bound the above by

Using a confidence parameter $\delta$, we pick $\beta=\gamma=\frac{\delta}{2}$, which ensures that the accuracy bound above holds with probability at least

for suitable positive constants $c,c^{\prime},c^{\prime\prime}$.

To get the theorem statement, recall that the analysis we performed pertains to data whose squared norm is bounded by $1$. By the reduction discussed at the beginning of the proof, we can apply it to data with squared norm at most $r$, by replacing $\lambda$ with $\lambda/r$, and $\eta$ with $\eta r$, leading to the condition

2 Proof of Theorem 2

The proof relies mainly on the techniques and lemmas of Section 6.1, used to prove Theorem 1. As done in Section 6.1, we will assume without loss of generality that $r=\max_{i}\|\mathbf{x}_{i}\|^{2}$ is at most $1$, and then transform the bound to a bound for general $r$ (see the discussion at the beginning of Subsection 6.1.2)

First, we extract the following result, which is essentially the first part of Lemma 11 (for $k=1$):

Let $A,\mathbf{w}_{t}$ be as defined in Algorithm 1, and suppose that $\eta\in\left[0,\frac{1}{23}\right]$. Then

for some positive numerical constants $c,c^{\prime}$.

The proof is based on martingale arguments, quite similar to the ones in Subsection 6.1.2 but with slight changes. First, we let

to simplify notation. We note that $b_{0}=1-\langle\mathbf{v}_{1},\mathbf{w}_{0}\rangle^{2}$ is assumed fixed, whereas $b_{1},b_{2},\ldots$ are random variables based on the sampling process. Lemma 11 tells us that if $\eta$ is sufficiently small, and $b_{t}\leq 1-\xi$ for some $\xi\in(0,1)$, then

for some numerical constants $c,c^{\prime}$.

Let $B$ be the event that $b_{t}\leq 1-\xi$ for all $t=0,1,\ldots,T$. Then for certain positive numerical constants $c_{1},c_{2},c_{3}$, if $\eta\leq c_{1}\lambda$, then

Using Eq. (21), we have for any $b_{t}$ satisfying event $B$ that

Taking expectation over $b_{t}$ (conditioned on $B$), we get that

We now turn to prove that the event $B$ assumed in Lemma 12 indeed holds with high probability:

The following holds for certain positive numerical constants $c_{1},c_{2},c_{3}$: If $\eta\leq c_{1}\lambda$, then for any $\beta\in(0,1)$, if

then it holds with probability at least $1-\beta$ that

To prove the lemma, we analyze the stochastic process $b_{1},b_{2},\ldots,b_{T}$, and use a concentration of measure argument. First, we collect the following facts:

$b_{0}\leq 1-\xi$: This directly follows from the assumption stated in the lemma.

$|b_{t+1}-b_{t}|$ is bounded by $c\eta$ for some constant $c$: Applying Lemma 9 for the case $k=1$, and assuming $\eta\leq 1/12$,

Armed with these facts, and using the maximal version of the Hoeffding-Azuma inequality , it follows that with probability at least $1-\beta$, it holds simultaneously for all $t=0,1,\ldots,T$ that

for some constants $c_{2},c_{3}$. If the expression is indeed less than $1-\xi$, then we get that $b_{t}\leq 1-\xi$ for all $t$, from which the lemma follows. ∎

Combining Lemma 16 and Lemma 17, and using Markov’s inequality, we get the following corollary:

Let confidence parameters $\beta,\gamma\in(0,1)$ be fixed. Then for some positive numerical constants $c_{1},c_{2},c_{3},c,c^{\prime}$, if $\eta\leq c_{1}\lambda$ and

then with probability at least $1-(\beta+\gamma)$, it holds that

We are now ready to prove our theorem. By Lemma 18, for any $\beta,\gamma\in\left(0,\frac{1}{2}\right)$ and any

we get with probability at least $1-(\beta+\gamma)$ that

Using the inequality $(1-(1/x))^{ax}\leq\exp(-a)$, which holds for any $x>1$ and any $a$, and taking $x=1/(c\eta\lambda\xi)$ and $a=3\log(1/\gamma)$, we can upper bound the above by

and since we assume $\gamma<\frac{1}{2}$, this is at most $\frac{1}{2}$. Overall, we got that with probability at least $1-\beta-\gamma$, $b_{T}\leq\frac{1}{2}$, and therefore $1-\langle\mathbf{v}_{1},\mathbf{w}_{T}\rangle^{2}\leq\frac{1}{2}$ as required.

It remains to show that the parameter choices in Eq. (23) can indeed be satisfied. First, we fix $\xi=\frac{1}{2}\zeta$ (where we recall that $0<\zeta\leq\langle\mathbf{v}_{1},\mathbf{w}_{0}\rangle^{2}$), which trivially ensures that $b_{0}=1-\langle\mathbf{v}_{1},\mathbf{w}_{0}\rangle^{2}$ is at most $1-2\xi$. Moreover, suppose we pick $\beta=\gamma$ in $(0,\exp(-1))$, and $\eta,T$ so that

where $c_{*},c^{\prime}_{*}$ are sufficiently small constants so that the bounds on $\eta,T$ in Eq. (23) are satisfied. This implies that the third bound in Eq. (23) is also satisfied, since by plugging in the values / bounds of $T$ and $\eta$, and using the assumptions $\gamma=\beta\leq\exp(-1)$ and $\xi\leq 1$, we have

which is less than $1-\xi$ if we pick $c_{*}$ sufficiently small compared to $c^{\prime}_{*}$.

To summarize, we get that for any $\gamma\in(0,\exp(-1))$, by picking $\eta$ as in Eq. (24), we have that after $T$ iterations (where $T$ is specified in Eq. (24)), with probability at least $1-2\gamma$, we get $\mathbf{w}_{T}$ such that $1-\langle\mathbf{v}_{1},\mathbf{w}_{T}\rangle\leq\frac{1}{2}$. Substituting $\delta=2\gamma$ and $\zeta=2\xi$, we get that if

(for some universal constant $c_{1}$), then with probability at least $1-\delta$, after

stochastic iterations, we get a satisfactory point $\mathbf{w}_{T}$.

As discussed at the beginning of the proof, this analysis is valid assuming $r=\max_{i}\|\mathbf{x}_{i}\|^{2}\leq 1$. By the reduction discussed at the beginning of Subsection 6.1.2, we can get an analysis for any $r$ by substituting $\lambda\rightarrow\lambda/r$ and $\eta\rightarrow\eta r$. This means that we should pick $\eta$ satisfying

3 Proof of Theorem 4

We first prove the following two auxiliary lemmas:

If $A$ is a symmetric matrix, then the gradient of the function $F(\mathbf{w})=-\frac{\mathbf{w}^{\top}A\mathbf{w}}{\|\mathbf{w}\|^{2}}$ at some $\mathbf{w}$ equals

where $B^{\bot}=B+B^{\top}$ (i.e., a matrix $B$ plus its transpose).