Streaming PCA: Matching Matrix Bernstein and Near-Optimal Finite Sample Guarantees for Oja's Algorithm

Prateek Jain, Chi Jin, Sham M. Kakade, Praneeth Netrapalli, Aaron Sidford

Introduction

Principal component analysis (PCA) is one of the most fundamental problems in machine learning, numerical linear algebra, and data analysis. It is commonly used for data compression, image processing, and visualization etc.

When we desire to perform PCA on large data sets, it may be the case that we cannot afford more than single pass over the data (or worse to even store the data in the first place) . To alleviate this issue, a popular line of research over the past several decades has been to consider streaming algorithms for PCA under the assumption that the data has reasonable statistical properties . There have been significant breakthroughs in getting near-optimal streaming PCA algorithms under fairly specialized models, e.g. spiked covariance .

This work considers one of the most natural variants of PCA, estimating the top eigenvector of a symmetric matrix, under a mild (and standard) set of assumptions under which concentration of measure applies (under the matrix Bernstein inequality). In particular, the setting is as follows:

∥Ai−Σ∥2≤M\left\|{\mathbf{A}_{i}-\mathbf{\Sigma}}\right\|_{2}\leq\mathcal{M} with probability 11, and

Let v1,...,vd\mathbf{v}_{1},...,\mathbf{v}_{d} denote the eigenvectors of Σ\mathbf{\Sigma} and λ1≥...≥λd\lambda_{1}\geq...\geq\lambda_{d} denote the corresponding eigenvalues. Our goal is to compute an ϵ\epsilon-approximation to v1\mathbf{v}_{1}, that is a unit vector w\mathbf{w} such that sin⁡2(w,v1)=△1−(w⊤v1)2≤ϵ\sin^{2}(\mathbf{w},\mathbf{v}_{1})\stackrel{{\scriptstyle\triangle}}{{=}}1-(\mathbf{w}^{\top}\mathbf{v}_{1})^{2}\leq\epsilon, in a single pass while minimizing space, time, and error (i.e. ϵ\epsilon). Note that sin⁡(w,v1)\sin(\mathbf{w},\mathbf{v}_{1}) denotes the sin⁡\sin of the angle between w\mathbf{w} and v1\mathbf{v}_{1}.

It is well known that to solve the Streaming PCA problem, one can simply compute the empirical covariance matrix 1n∑i∈[n]Ai\frac{1}{n}\sum_{i\in[n]}\mathbf{A}_{i} and compute the right singular vector of this matrix. Here, matrix Bernstein inequality and Wedin’s theorem implies the following standard sample complexity bound for the Streaming PCA problem:

Under the assumptions of Definition 1, the top right singular vector v^\widehat{\mathbf{v}} of Σ^=1n∑i∈[n]Ai\widehat{\mathbf{\Sigma}}=\frac{1}{n}\sum_{i\in[n]}\mathbf{A}_{i} is an ϵ\epsilon-approximation to the top eigenvector v1\mathbf{v}_{1} of Σ\mathbf{\Sigma} with probability 1−δ1-\delta, where

Theorem 1.1 is essentially the previous best sample complexity known for estimating the top eigenvector In recent work in it was shown that the log(d/δ)log(d/\delta) factor in the first term could be removed asymptotically for small enough ϵ\epsilon if only constant success probability is required.. Unfortunately, the above is purely a statistical claim, and, algorithmically, there are least two concerns. First, computing the empirical covariance matrix Σ^=1n∑i∈[n]Ai\widehat{\mathbf{\Sigma}}=\frac{1}{n}\sum_{i\in[n]}\mathbf{A}_{i} naively requires O(d2)O(d^{2}) time and space, and second, computing the top eigenvector of the empirical covariance matrix in general may require super linear time. While there have been many attempts to produce streaming algorithms that use only O(d)O(d) space to solve the streaming PCA problem, to our knowledge, all previous methods either lose a multiplicative factor of either λ1λ1−λ2\frac{\lambda_{1}}{\lambda_{1}-\lambda_{2}} or dd in the analysis in order to achieve constant accuracy when applied in our setting.

In an attempt to overcome this limitation and improve the guarantees for solving the streaming PCA problem, this work seeks to address the following question:

Can we match the sample complexity of matrix Bernstein + Wedin’s theorem with an algorithm that uses O(d)O(d) space only and takes a single linear-time pass over the input?

This work answers this question in the affirmative, showing that one can succeed with constant probability matching the sample complexity of Theorem 1.1 up to logarithmic terms and small additive factors. Interestingly, this is achieved by providing a novel analysis of the classical Oja’s algorithm, which is perhaps, the most popular algorithm for Streaming PCA.

This work shows that for proper choice of learning rates ηi\eta_{i}, Oja’s algorithm in fact can improve the best known results for streaming PCA and answer our question in the affirmative. In particular, we have that:

Let the assumptions of Definition 1 hold. Suppose the step size sequence for Algorithm 1 is chosen to be ηi=log⁡d(λ1−λ2)(β+i)\eta_{i}=\frac{\log d}{(\lambda_{1}-\lambda_{2})(\beta+i)}, where

Then the output wn\mathbf{w_{n}} of Algorithm 1 is an ϵ\epsilon-approximation to the top eigenvector v1\mathbf{v}_{1} of Σ\mathbf{\Sigma} satisfying

with probability greater than 3/43/4. Here CC is an absolute numerical constant.

The error above should be interpreted as being the sum of a O(1n)\mathcal{O}\left(\frac{1}{n}\right) higher order term and another O((2β/n)2log⁡d)\mathcal{O}\left((2\beta/n)^{2\log d}\right) lower order term which is at most o(1nlog⁡d)o(\frac{1}{n^{\log d}}) (once n>4β2n>4\beta^{2}). In particular, this result shows that, up to an additive lower order term, one can match Theorem 1.1 with an asymptotic error of O(Vlog⁡d(λ1−λ2)2n)\mathcal{O}\left(\frac{\mathcal{V}\log d}{(\lambda_{1}-\lambda_{2})^{2}n}\right) with constant probability. The lower order term has β\beta which is the max⁡\max of three parts: Mlog⁡d(λ1−λ2)\frac{\mathcal{M}\log d}{(\lambda_{1}-\lambda_{2})}, Vlog⁡2d(λ1−λ2)2\frac{\mathcal{V}\log^{2}d}{(\lambda_{1}-\lambda_{2})^{2}} and λ12log⁡2d(λ1−λ2)2\frac{\lambda_{1}^{2}\log^{2}d}{(\lambda_{1}-\lambda_{2})^{2}}. The first part, depending on M\mathcal{M}, is exactly the same as what appears in Theorem 1.1. The second one, depending on V\mathcal{V} has an additional log⁡d\log d factor over the first order term and is irrelevant once, say n>10βn>10\beta. Notably, the third part, depending on λ12\lambda_{1}^{2}, does not appear in Theorem 1.1; it arises here entirely due to computational reasons: the setting allows only a single linear-time pass over the matrices, while Theorem 1.1 makes no such assumption. For instance, consider the case V=0\mathcal{V}=0 which means A1=Σ\mathbf{A}_{1}=\mathbf{\Sigma}. Matrix Bernstein tells us that one sample is sufficient to compute v1\mathbf{v}_{1}. However, it is not evident how to compute it using a single pass over A1\mathbf{A}_{1}. Note however, that the rate at which the lower order terms, i.e. o(1nlog⁡d)o\left(\frac{1}{n^{\log d}}\right), decrease is much better than O(1/n2)\mathcal{O}\left(1/n^{2}\right) guaranteed by Theorem 1.1.

In fact, this result also improves the asymptotic error rate obtained by Theorem 1.1. In particular, the following result shows that Oja’s algorithm gets an asymptotic rate of O(V(λ1−λ2)2n)\mathcal{O}\left(\frac{\mathcal{V}}{(\lambda_{1}-\lambda_{2})^{2}n}\right) which is better than that of matrix Bernstein by a factor of O(log⁡d)\mathcal{O}\left(\log d\right).A similar asymptotic result was recently obtained by. However, their result requires an initial vector that is constant close to v1\mathbf{v}_{1}, which itself is a difficult problem.

Let the assumptions of Definition 1 hold. Suppose the step size sequence for Algorithm 1 is chosen to be ηi=6(λ1−λ2)(β+i)\eta_{i}=\frac{6}{(\lambda_{1}-\lambda_{2})(\beta+i)}, where

Suppose n>β1.2d0.1n>\beta^{1.2}d^{0.1}. Then the output wn\mathbf{w_{n}} of Algorithm 1 is an ϵ\epsilon-approximation to the top eigenvector v1\mathbf{v}_{1} of Σ\mathbf{\Sigma} satisfying

with probability greater than 3/43/4. Here CC is an absolute numerical constant.

Note that Theorems 1.2 and 1.3 guarantee success probability of 3/43/4. One way to boost the probability to 1−δ1-\delta, for some δ>0\delta>0, is to run O(log⁡1/δ)\mathcal{O}\left(\log 1/\delta\right) copies of the algorithm, each with 3/43/4 success probability and then output the geometric median of the solutions, which can be done in nearly linear time. The detailes are omitted here.

Beyond the improved sample complexities we believe our analysis sheds light on the type of step sizes for which Oja’s algorithm converges quickly and therefore illuminates how to efficiently perform streaming PCA. We note that we have essentially assumed an oracle which sets the step size sequence, and an important question is how to set the step size in a robust and data data driven manner. Moreover, we believe that our analysis is fairly general and hope that it may be extended to make progress on analyzing the many variants of PCA that occur in both theory and in practice.

Here we compare our sample complexity bounds with existing analyses of various methods. Recall that the error of the estimate w\mathbf{w} is sin⁡2(w,v1)=1−(w⊤v1)2\sin^{2}(\mathbf{w},\mathbf{v}_{1})=1-(\mathbf{w}^{\top}\mathbf{v}_{1})^{2}.

We consider three popular methods used for computing v1\mathbf{v}_{1}. The first one is the batch method which computes largest eigenvector of empirical covariance and uses Wedin’s theorem with matrix Bernstein inequality (cf. Theorem 1.1). The second method is Alecton, which is very similar to Oja’s algorithm . Finally, consider a block-power method (BPM) which divides samples into different blocks and applies power iteration to the empirical estimate from each block. See Table 1 for the comparison.

We stress that some of the results we compare to make different assumptions than Definition 1. The bounds stated for them are our best attempt to adapt their bounds in the setting of Definition 1 (which is quite standard). The next paragraph provides a simple example, which demonstrates the improvement in our result as compared to existing work.

2 Additional Related Work

Existing results for computing largest eigenvector of a data covariance matrix using streaming samples can be divided into three broad settings: a) stochastic data, b) arbitrary sequence of data, c) regret bounds for arbitrary sequence of data.

Stochastic data: Here, the data is assumed to be sampled i.i.d. from a fixed distribution. The analysis of Oja’s algorithm as well as those of block power method and Alecton mentioned earlier are in this setting. also obtained a result in the restricted spiked covariance model. provides an analysis of a modification of Oja’s algorithm but with an extra O(d5)O(d^{5}) multiplicative factor compared to ours. provides an algorithm based on shift and invert framework that obtains the same asymptotic error as ours. However, their algorithm requires warm start with a vector that is already constant close to the top eigenvector, which itself is a hard problem.

Arbitrary data: In this setting, each row of the data matrix is provided in an arbitrary order. Most of the existing methods here first compute a sketch of the matrix and use that to compute an estimate of the top eigenvector . However, a direct application of such techniques to the stochastic setting leads to sample complexity bounds which are larger by a multiplicative factor of O(d)O(d) (ignoring other factors like variance etc). Finally, also provide methods for eigenvector computation, but they require multiple passes over the data and hence do not apply to the streaming setting.

Regret bounds: Here, at each step the algorithm has to output an estimate w\mathbf{w} of v1\mathbf{v}_{1} for which we get reward of wTAiw\mathbf{w}^{T}\mathbf{A}_{i}\mathbf{w} and the goal is to minimize the regret w.r.t. v1\mathbf{v}_{1}. The algorithms in this regime are mostly based on online convex optimization and applying them in our setting would again result in a loss of multiplicative O(d)O(d). Moreover, typical algorithms in this setting are not memory efficient .

3 Notation

4 Paper Organization

The rest of this paper is organized as follows. Section 2 introduces basic mathematical facts used throughout the paper and also provides a proof of the error bound of the standard batch method (Theorem 1.1). Section 3 provides an overview of our approach to analyzing Oja’s algorithm and provides the main technical result of the paper. This technical result is used in Section 4 to prove the running time for Oja’s algorithm and to justify the choice of step size. Section 5 presents the proof of the main technical result. Section 6 concludes and mentions a few interesting future directions.

Preliminaries

The following basic inequalities regarding power series, the exponential, and PSD matrices are used throughout. The facts are summarized here:

1+x≥exp⁡(x−x2)1+x\geq\exp\left(x-x^{2}\right) for all x≥0x\geq 0

11+x≤∑i=1∞1(x+i)2≤1x\frac{1}{1+x}\leq\sum_{i=1}^{\infty}\frac{1}{(x+i)^{2}}\leq\frac{1}{x}

⟨A,B⟩≤⟨A,C⟩\left\langle\mathbf{A},\mathbf{B}\right\rangle\leq\left\langle\mathbf{A},\mathbf{C}\right\rangle for PSD matrices A,B,C\mathbf{A},\mathbf{B},\mathbf{C} with B⪯C\mathbf{B}\preceq\mathbf{C}

The first inequality follows from the Taylor expansion of exp⁡(x)\exp(x). The second comes from 1+0=exp⁡(0−02)1+0=\exp(0-0^{2}) and ddx(1+x)≤ddxexp⁡(x−x2)\frac{d}{dx}(1+x)\leq\frac{d}{dx}\exp(x-x^{2}) for x≥0x\geq 0. The third follows by considering upper and lower Riemann sums of ∫y=1∞1/(x+y)\int_{y=1}^{\infty}1/(x+y). The fourth from the fact that since A\mathbf{A} is PSD there is a matrix D\mathbf{D} with D⊤D=A\mathbf{D}^{\top}\mathbf{D}=\mathbf{A} and therefore

The final follows from Cauchy Schwarz and Young’s inequality, i.e. x⋅y≤12(x2+y2)x\cdot y\leq\frac{1}{2}(x^{2}+y^{2}) as

The following is a matrix Bernstein based proof of the error bound of the batch method.

Using Theorem 1.4 of , we have (w.p. ≥1−δ\geq 1-\delta):

Let v^\widehat{\mathbf{v}} be the top eigenvector of Σ^=1n∑i=1nAi\widehat{\mathbf{\Sigma}}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{A}_{i}. Using Wedin’s theorem , implies:

Theorem now follows by combining (1) and (2). ∎

Approach

Let us now describe the approach to analyze Oja’s algorithm. We provide our main theorem regarding the convergence rate of Oja’s algorithm and discuss how it is proved. The details of the proof are deferred to Section 5 and the use of the theorem to choose step sizes is in Section 4.

One of the primary difficulties in analyzing Oja’s algorithm, or more broadly any algorithm for streaming PCA, is choosing a subtle potential function to analyze the method. If we try to analyze the progress of Oja’s algorithm in every iteration ii, by measuring the quality of wi\mathbf{w}_{i}, we run the risk that during the first few iterations of Oja’s algorithm a step may actually yield a wi+1\mathbf{w}_{i+1} that is orthogonal to vi\mathbf{v}_{i}. If this happens, even in the typical best case, where all future samples are Σ\mathbf{\Sigma} itself, we would still fail to converge. In short, if we do not account for the randomness of w0\mathbf{w}_{0} in our potential function then it is difficult to show that a rapidly convergent algorithm does not catastrophically fail.

Rather than analyzing the convergence of wi\mathbf{w}_{i} directly we instead analyze the convergence of Oja’s algorithm as an operator on w0\mathbf{w}_{0}. Oja’s algorithm simply considers the matrix

and outputs the normalized result of applying this matrix, Bn\mathbf{B}_{n}, to the random initial vector, i.e.

Rather than analyze the improvement of wn+1\mathbf{w}_{n+1} over wn\mathbf{w}_{n} we analyze Bn+1\mathbf{B}_{n+1}’s improvement over Bn\mathbf{B}_{n}.

Another interpretation of (3) and (4) is that Oja’s algorithm simply approximates vn\mathbf{v}_{n} by performing 1 step of the power method on the matrix Bn\mathbf{B}_{n}. Fortunately, analyzing when 1 step of the power method succeeds is fairly straightforward as we show below:

As w\mathbf{w} is distributed uniformly over the sphere, we have: w=g/∥g∥2\mathbf{w}=\mathbf{g}/\|\mathbf{g}\|_{2} where g∼N(0,I)\mathbf{g}\sim N(0,I). Consequently, with probability at least 1−δ1-\delta

Let δ>0\delta>0 and step sizes ηi≤14⋅max⁡{M,λ1}\eta_{i}\leq\frac{1}{4\cdot\max\{M,\lambda_{1}\}}. The output wn\mathbf{w_{n}} of Algorithm 1 is an ϵ\epsilon-approximation to v1\mathbf{v}_{1} with probability at least 1−δ1-\delta where

where Q=△δ2Clog⁡(1/δ)(1−1δexp⁡(18V‾∑i=1nηi2)−1)Q\stackrel{{\scriptstyle\triangle}}{{=}}\frac{\delta^{2}}{C\log(1/\delta)}\left(1-\frac{1}{\sqrt{\delta}}\sqrt{\exp\left(18\overline{\mathcal{V}}\sum_{i=1}^{n}\eta_{i}^{2}\right)-1}\right), V‾=△V+λ12\overline{\mathcal{V}}\stackrel{{\scriptstyle\triangle}}{{=}}\mathcal{V}+\lambda_{1}^{2}, and CC is an absolute constant.

Theorem 3.1 is proved in Section 5. Theorem 3.1 serves as the basis for our results regarding Oja’s algorithm. In the next section we show how to use this theorem to choose step sizes and achieve the main results of this paper.

Main Results

Theorem 3.1, from the previous section, leads to our main results, provided here. The theorem and proof are below and essentially consist of choosing appropriate parameters to efficiently apply Theorem 3.1. Once we have this theorem, Theorems 1.2 and 1.3 follow by choosing α=log⁡d\alpha=\log d and α=6\alpha=6 respectively.

Fix any δ>0\delta>0 and suppose the step sizes are set to ηt=α(λ1−λ2)(β+t)\eta_{t}=\frac{\alpha}{(\lambda_{1}-\lambda_{2})(\beta+t)} for α>12\alpha>\frac{1}{2} and

Suppose the number of samples n>βn>\beta. Then the output wn\mathbf{w_{n}} of Algorithm 1 satisfies:

with probability at least 1−δ1-\delta. Here CC is an absolute numerical constant.

where Q=△δ2Clog⁡(1/δ)(1−1δexp⁡(18V‾∑i=1nηi2)−1)Q\stackrel{{\scriptstyle\triangle}}{{=}}\frac{\delta^{2}}{C\log(1/\delta)}\left(1-\frac{1}{\sqrt{\delta}}\sqrt{\exp\left(18\overline{\mathcal{V}}\sum_{i=1}^{n}\eta_{i}^{2}\right)-1}\right). Since ηi=α(λ1−λ2)(β+i)\eta_{i}=\frac{\alpha}{(\lambda_{1}-\lambda_{2})(\beta+i)}, we have ∑i∈[n]ηi2≤α2(λ1−λ2)2β\sum_{i\in[n]}\eta_{i}^{2}\leq\frac{\alpha^{2}}{(\lambda_{1}-\lambda_{2})^{2}\beta} and by our assumption that V‾α2(λ1−λ2)2β≤118log⁡(1+δ100)\frac{\overline{\mathcal{V}}\alpha^{2}}{(\lambda_{1}-\lambda_{2})^{2}\beta}\leq\frac{1}{18}\log\left(1+\frac{\delta}{100}\right), we have:

Moreover, since ∑i∈[n]ηi≥αλ1−λ2log⁡(1+n/β)\sum_{i\in[n]}\eta_{i}\geq\frac{\alpha}{\lambda_{1}-\lambda_{2}}\log\left(1+n/\beta\right), we have

Note that ∑j=i+1nηj≤αλ1−λ2log⁡n+β+1i+β+1\sum_{j=i+1}^{n}\eta_{j}\leq\frac{\alpha}{\lambda_{1}-\lambda_{2}}\log\frac{n+\beta+1}{i+\beta+1}. Moreover, as α>1/2\alpha>1/2, we have:

Substituting (6), (7) and (8) into (5) proves the theorem. ∎

Bounding the Convergence of Oja’s Algorithm

In this section, we present a detailed proof of Theorem 3.1. The proof follows the approach outlined in Section 3 and uses the notation of that section, i.e.

We let Bn=△(I+ηnAn)⋯(I+η1A1)\mathbf{B}_{n}\stackrel{{\scriptstyle\triangle}}{{=}}\left(\mathbf{I}+\eta_{n}\mathbf{A}_{n}\right)\cdots\left(\mathbf{I}+\eta_{1}\mathbf{A}_{1}\right) with B0=△I\mathbf{B}_{0}\stackrel{{\scriptstyle\triangle}}{{=}}\mathbf{I}

We let V‾=△V+λ12\overline{\mathcal{V}}\stackrel{{\scriptstyle\triangle}}{{=}}\mathcal{V}+\lambda_{1}^{2}

For all t≥0t\geq 0 and ηi≥0\eta_{i}\geq 0 we have

The result follows by using induction along with α0=1\alpha_{0}=1 and 1+x≤ex1+x\leq e^{x}. ∎

For all t≥0t\geq 0 and ηi≤1λ1\eta_{i}\leq\frac{1}{\lambda_{1}} the following holds

where ζ1\zeta_{1} follows from the fact that V⊥\mathbf{V}_{\perp} is orthogonal to v1\mathbf{v}_{1} and ζ2\zeta_{2} follows from defintion of V\mathcal{V}.

Plugging the above into (10), we get for all t≥1t\geq 1,

where the last inequality follows from 1+x≤ex1+x\leq e^{x} and using Lemma 5.1.

Recursing the above inequality, we obtain

Since B0=I\mathbf{B}_{0}=\mathbf{I} we see that α0=d−1≤d\alpha_{0}=d-1\leq d. Using that ηi≤1λ1≤1λ2\eta_{i}\leq\frac{1}{\lambda_{1}}\leq\frac{1}{\lambda_{2}} completes the proof. ∎

For all t≥0t\geq 0 and ηi≥0\eta_{i}\geq 0 we have

Consequently βt≥(1+2ηtλ1)βt−1\beta_{t}\geq(1+2\eta_{t}\lambda_{1})\beta_{t-1}. Furthermore, B0=I\mathbf{B}_{0}=\mathbf{I} and hence β0=∥v1∥22=1\beta_{0}=\left\|{\mathbf{v}_{1}}\right\|_{2}^{2}=1. Proceeding by induction and using that 1+x≥exp⁡(x−x2)1+x\geq\exp(x-x^{2}) for all x≥0x\geq 0 finishes the proof. ∎

For t≥0t\geq 0 suppose that ηi≤14⋅max⁡{λ1,M}\eta_{i}\leq\frac{1}{4\cdot\max\{\lambda_{1},M\}} for all i∈[t]i\in[t] then.

where Gt−1=△Wt,t−1⊤v1v1⊤Wt,t−1{\mathbf{G}_{t-1}}\stackrel{{\scriptstyle\triangle}}{{=}}{\mathbf{W}_{t,t-1}^{\top}\mathbf{v}_{1}\mathbf{v}_{1}^{\top}\mathbf{W}_{t,t-1}}. In order to bound the above quantity, we first bound the above expression for an arbitrary Gt−1≡G{\mathbf{G}_{t-1}}\equiv\mathbf{G}. We then take an expectation over only A1\mathbf{A}_{1} and then finally take an expectation over Gt−1{\mathbf{G}_{t-1}}. That is, for an arbitrary fixed symmetric matrix G\mathbf{G}, we have:

We now bound the various terms above as follows. Each of the second order terms can be bounded using Lemma 2.1 as follows:

The third order terms can be bounded as follows:

where we used the assumption that ∥A1∥2≤∥A1−Σ∥2+∥Σ∥2≤M+λ1\left\|{\mathbf{A}_{1}}\right\|_{2}\leq\left\|{\mathbf{A}_{1}-\mathbf{\Sigma}}\right\|_{2}+\left\|{\mathbf{\Sigma}}\right\|_{2}\leq\mathcal{M}+\lambda_{1} with probability 11. Finally the fourth order term can be bounded as

Plugging (13), (14) and (15) into (12) tells us that

where in the last line we used that ηi≤14max⁡{M,λ1}\eta_{i}\leq\frac{1}{4\max\{\mathcal{M},\lambda_{1}\}} and that 1+x≤exp⁡(x)1+x\leq\exp(x)

Using the value G=Gt−1=Wt,t−1⊤v1v1⊤Wt,t−1\mathbf{G}={\mathbf{G}_{t-1}}={\mathbf{W}_{t,t-1}^{\top}\mathbf{v}_{1}\mathbf{v}_{1}^{\top}\mathbf{W}_{t,t-1}} and plugging the above into (11), we have

We now have everything to prove Theorem 3.1.

First, using Chebyshev’s inequality, we have:

So with probability greater than 1−δ1-\delta, the following holds:

where ζ1\zeta_{1} follows from Lemma 5.3 and 5.4.

Furthermore, using Lemma 5.2 and Markov’s inequality, we have with probability at least 1−δ1-\delta,

Consequently with probability at least 1−2δ1-2\delta both (LABEL:eqn:main1) and (17) hold and therefore the result follows by Lemma 3.1 and choosing a δ\delta that is smaller by a constant. ∎

Conclusion and Future Work

This work presented a finite sample complexity and asymptotic convergence rates for the classic Oja’s algorithm for top-11 component streaming PCA that match well known matrix concentration and perturbation results for computing the top eigenvector. In fact, asymptotically our bound improves upon standard matrix Bernstein bounds by a factor of O(log⁡d)\mathcal{O}\left(\log d\right). Our results are tighter than existing streaming PCA results by a factor of either O(d)\mathcal{O}\left(d\right) or O(1/gap)\mathcal{O}\left(1/\textrm{gap}\right).

Our analysis relied on a novel view of the algorithm and is technically fairly simple. We hope that our analysis opens a way to make progress on the many variants of PCA that occur in both theory and practice. In particular, we believe the following directions should be of wide interest:

Multiple components: Currently, our result holds only for estimating the top eigenvector of Σ\mathbf{\Sigma}. Extension of our technique to compute top-kk eigenvectors is an important future direction.

Rayleigh quotient: Another standard metric to measure optimality of wn\mathbf{w_{n}} is Rayleigh quotient: wn⊤Σwn\mathbf{w_{n}}^{\top}\mathbf{\Sigma}\mathbf{w_{n}}. Converting our bounds on sin⁡2(wn,v1)\sin^{2}(\mathbf{w_{n}},\mathbf{v}_{1}) to Rayleigh quotient loses a multiplicative factor of O(1/gap)\mathcal{O}\left(1/\textrm{gap}\right) compared to the optimal rate. A direct analysis that does not lose this factor is an interesting open problem. Results on Rayleigh quotient may also help in obtaining sample complexity guarantees that are independent of eigenvalue gap.

High Probability: This work focused on obtaining tight bounds on the error. However, the dependence of our results on success probability is quite suboptimal. One way to fix this is to run many copies of the algorithm, each with say 3/43/4 success probability and then output the geometric median of the solutions, which can be done in nearly linear time. However, we conjecture that a tighter analysis using our techniques might directly lead to improved dependency on success probability and possibly help solve some of the other problems mentioned above.

Acknowledgements

Sham Kakade acknowledges funding from the Washington Research Foundation for innovation in Data-intensive Discovery.

References