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:
with probability , and
Let denote the eigenvectors of and denote the corresponding eigenvalues. Our goal is to compute an -approximation to , that is a unit vector such that , in a single pass while minimizing space, time, and error (i.e. ). Note that denotes the of the angle between and .
It is well known that to solve the Streaming PCA problem, one can simply compute the empirical covariance matrix 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 of is an -approximation to the top eigenvector of with probability , 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 factor in the first term could be removed asymptotically for small enough 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 naively requires 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 space to solve the streaming PCA problem, to our knowledge, all previous methods either lose a multiplicative factor of either or 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 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 , 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 , where
Then the output of Algorithm 1 is an -approximation to the top eigenvector of satisfying
with probability greater than . Here is an absolute numerical constant.
The error above should be interpreted as being the sum of a higher order term and another lower order term which is at most (once ). In particular, this result shows that, up to an additive lower order term, one can match Theorem 1.1 with an asymptotic error of with constant probability. The lower order term has which is the of three parts: , and . The first part, depending on , is exactly the same as what appears in Theorem 1.1. The second one, depending on has an additional factor over the first order term and is irrelevant once, say . Notably, the third part, depending on , 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 which means . Matrix Bernstein tells us that one sample is sufficient to compute . However, it is not evident how to compute it using a single pass over . Note however, that the rate at which the lower order terms, i.e. , decrease is much better than 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 which is better than that of matrix Bernstein by a factor of .A similar asymptotic result was recently obtained by. However, their result requires an initial vector that is constant close to , 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 , where
Suppose . Then the output of Algorithm 1 is an -approximation to the top eigenvector of satisfying
with probability greater than . Here is an absolute numerical constant.
Note that Theorems 1.2 and 1.3 guarantee success probability of . One way to boost the probability to , for some , is to run copies of the algorithm, each with 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 is .
We consider three popular methods used for computing . 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 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 (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 of for which we get reward of and the goal is to minimize the regret w.r.t. . 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 . 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:
for all
for PSD matrices with
The first inequality follows from the Taylor expansion of . The second comes from and for . The third follows by considering upper and lower Riemann sums of . The fourth from the fact that since is PSD there is a matrix with and therefore
The final follows from Cauchy Schwarz and Young’s inequality, i.e. 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. ):
Let be the top eigenvector of . 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 , by measuring the quality of , we run the risk that during the first few iterations of Oja’s algorithm a step may actually yield a that is orthogonal to . If this happens, even in the typical best case, where all future samples are itself, we would still fail to converge. In short, if we do not account for the randomness of 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 directly we instead analyze the convergence of Oja’s algorithm as an operator on . Oja’s algorithm simply considers the matrix
and outputs the normalized result of applying this matrix, , to the random initial vector, i.e.
Rather than analyze the improvement of over we analyze ’s improvement over .
Another interpretation of (3) and (4) is that Oja’s algorithm simply approximates by performing 1 step of the power method on the matrix . Fortunately, analyzing when 1 step of the power method succeeds is fairly straightforward as we show below:
As is distributed uniformly over the sphere, we have: where . Consequently, with probability at least
Let and step sizes . The output of Algorithm 1 is an -approximation to with probability at least where
where , , and 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 and respectively.
Fix any and suppose the step sizes are set to for and
Suppose the number of samples . Then the output of Algorithm 1 satisfies:
with probability at least . Here is an absolute numerical constant.
where . Since , we have and by our assumption that , we have:
Moreover, since , we have
Note that . Moreover, as , 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 with
We let
For all and we have
The result follows by using induction along with and . ∎
For all and the following holds
where follows from the fact that is orthogonal to and follows from defintion of .
Plugging the above into (10), we get for all ,
where the last inequality follows from and using Lemma 5.1.
Recursing the above inequality, we obtain
Since we see that . Using that completes the proof. ∎
For all and we have
Consequently . Furthermore, and hence . Proceeding by induction and using that for all finishes the proof. ∎
For suppose that for all then.
where . In order to bound the above quantity, we first bound the above expression for an arbitrary . We then take an expectation over only and then finally take an expectation over . That is, for an arbitrary fixed symmetric matrix , 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 with probability . 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 and that
Using the value 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 , the following holds:
where follows from Lemma 5.3 and 5.4.
Furthermore, using Lemma 5.2 and Markov’s inequality, we have with probability at least ,
Consequently with probability at least both (LABEL:eqn:main1) and (17) hold and therefore the result follows by Lemma 3.1 and choosing a 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- 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 . Our results are tighter than existing streaming PCA results by a factor of either or .
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 . Extension of our technique to compute top- eigenvectors is an important future direction.
Rayleigh quotient: Another standard metric to measure optimality of is Rayleigh quotient: . Converting our bounds on to Rayleigh quotient loses a multiplicative factor of 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 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.