The Fast Convergence of Incremental PCA

Introduction

Here $\gamma_{n}$ is a “learning rate” that is typically proportional to $1/n$.

$\sum_{n}\gamma_{n}=\infty$ while $\sum_{n}\gamma_{n}^{2}<\infty$.

If $\lambda_{1},\lambda_{2}$ denote the top two eigenvalues of $A$, then $\lambda_{1}>\lambda_{2}$.

There are also other incremental estimators for which convergence has not been established; see, for instance, and .

In this paper, we analyze the rate of convergence of the Krasulina and Oja estimators. They can be treated in a common framework, as stochastic approximation algorithms for maximizing the Rayleigh quotient

The maximum value of this function is $\lambda_{1}$, and is achieved at $v^{*}$ (or any nonzero multiple thereof). The gradient is

Recently, there has been a lot of work on rates of convergence for stochastic gradient descent (for instance, ), but this has typically been limited to convex cost functions. These results do not apply to the non-convex Rayleigh quotient, except at the very end, when the system is near convergence. Most of our analysis focuses on the buildup to this finale.

We measure the quality of the solution $V_{n}$ at time $n$ using the potential function

Set starting time. Set the clock to time $n_{o}$.

Update step. Perform either the Krasulina or Oja update, with $\gamma_{n}=c/n$.

The first step is similar to using a learning rate of the form $\gamma_{n}=c/(n+n_{o})$, as is often done in stochastic gradient descent implementations . We have adopted it because the initial sequence of updates is highly noisy: during this phase $V_{n}$ moves around wildly, and cannot be shown to make progress. It becomes better behaved when the step size $\gamma_{n}$ becomes smaller, that is to say when $n$ gets larger than some suitable $n_{o}$. By setting the start time to $n_{o}$, we can simply fast-forward the analysis to this moment.

2 Initialization

One possible initialization is to set $V_{n_{o}}$ to the first data point that arrives, or to the average of a few data points. This seems sensible enough, but can fail dramatically in some situations.

Here is an example. Suppose $X$ can take on just $2d$ possible values: $\pm e_{1},\pm\sigma e_{2},\ldots,\pm\sigma e_{d}$, where the $e_{i}$ are coordinate directions and $0<\sigma<1$ is a small constant. Suppose further that the distribution of $X$ is specified by a single positive number $p<1$:

Then $X$ has mean zero and covariance $\mbox{diag}(p,\sigma^{2}(1-p)/(d-1),\ldots,\sigma^{2}(1-p)/(d-1))$. We will assume that $p$ and $\sigma$ are chosen so that $p>\sigma^{2}(1-p)/(d-1)$; in our notation, the top eigenvalues are then $\lambda_{1}=p$ and $\lambda_{2}=\sigma^{2}(1-p)/(d-1)$, and the target vector is $v^{*}=e_{1}$.

If $V_{n}$ is ever orthogonal to some $e_{i}$, it will remain so forever. This is because both the Krasulina and Oja updates have the following properties:

If $V_{n_{o}}$ is initialized to a random data point, then with probability $1-p$, it will be assigned to some $e_{i}$ with $i>1$, and will converge to a multiple of that same $e_{i}$ rather than to $e_{1}$. Likewise, if it is initialized to the average of $\leq 1/p$ data points, then with constant probability it will be orthogonal to $e_{1}$ and remain so always.

Setting $V_{n_{o}}$ to a random unit vector avoids this problem. However, there are doubtless cases, for instance when the data has intrinsic dimension $\ll d$, in which a better initializer is possible.

3 The setting of the learning rate

In order to get a sense of what rates of convergence we might expect, let’s return to the example of a random vector $X$ with $2d$ possible values. In the Oja update $V_{n}=V_{n-1}+\gamma_{n}X_{n}X_{n}^{T}V_{n-1}$, we can ignore normalization if we are merely interested in the progress of the potential function $\Psi_{n}$. Since the $X_{n}$ correspond to coordinate directions, each update changes just one coordinate of $V$:

Recall that we initialize $V_{n_{o}}$ to a random vector from the unit sphere. For simplicity, let’s just suppose that $n_{o}=0$ and that this initial value is the all-ones vector (again, we don’t have to worry about normalization). On each iteration the first coordinate is updated with probability exactly $p=\lambda_{1}$, and thus

since $\gamma_{n}=c/n$. Likewise, for $i>1$,

If all goes according to expectation, then at time $n$,

(This is all very rough, but can be made precise by obtaining concentration bounds for $\ln V_{n,i}$.) From this, we can see that it is not possible to achieve a $O(1/n)$ rate unless $c\geq 1/(2(\lambda_{1}-\lambda_{2}))$. Therefore, we will assume this when stating our final results, although most of our analysis is in terms of general $\gamma_{n}$. An interesting practical question, to which we do not have an answer, is how one would empirically set $c$ without prior knowledge of the eigenvalue gap.

4 Nested sample spaces

For $n\geq n_{o}$, let $\mathcal{F}_{n}$ denote the sigma-field of all outcomes up to and including time $n$: $\mathcal{F}_{n}=\sigma(V_{n_{o}},X_{n_{o}+1},\ldots,X_{n})$. We start by showing that

To deal with this, we divide the analysis into epochs: the first takes $\Psi_{n}$ from $1-1/d$ to $1-2/d$, the second from $1-2/d$ to $1-4/d$, and so on until $\Psi_{n}$ finally drops below $1/2$. We use martingale large deviation bounds to bound the length of each epoch, and also to argue that $\Psi_{n}$ does not regress. In particular, we establish a sequence of times $n_{j}$ such that (with high probability)

The analysis of each epoch uses martingale arguments, but at the same time, assumes that $\Psi_{n}$ remains bounded above. Combining the two requires a careful specification of the sample space at each step. Let $\Omega$ denote the sample space of all realizations $(v_{n_{o}},x_{n_{o}+1},x_{n_{o}+2},\ldots)$, and $P$ the probability distribution on these sequences. For any $\delta>0$, we define a nested sequence of spaces $\Omega\supset\Omega_{n_{o}}^{\prime}\supset\Omega_{n_{o}+1}^{\prime}\supset\cdots$ such that each $\Omega_{n}^{\prime}$ is $\mathcal{F}_{n-1}$-measurable, has probability $P(\Omega_{n}^{\prime})\geq 1-\delta$, and moreover consists exclusively of realizations $\omega\in\Omega$ that satisfy the constraints (1) up to and including time $n-1$. We can then build martingale arguments by restricting attention to $\Omega_{n}^{\prime}$ when computing the conditional expectations of quantities at time $n$.

5 Main result

There is a constant $B$ such that $\|X_{n}\|^{2}\leq B$.

The eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{d}$ of $A$ satisfy $\lambda_{1}>\lambda_{2}$.

The step sizes are of the form $\gamma_{n}=c/n$.

Under these conditions, we get the following rate of convergence for the Krasulina update.

There are absolute constants $A_{o},A_{1}>0$ and $1<a<4$ for which the following holds. Pick any $0<\delta<1$, and any $c_{o}>2$. Set the step sizes to $\gamma_{n}=c/n$, where $c=c_{o}/(2(\lambda_{1}-\lambda_{2}))$, and set the starting time to $n_{o}\geq(A_{o}B^{2}c^{2}d^{2}/\delta^{4})\ln(1/\delta)$. Then there is a nested sequence of subsets of the sample space $\Omega\supset\Omega_{n_{o}}^{\prime}\supset\Omega_{n_{o}+1}^{\prime}\supset\cdots$ such that for any $n\geq n_{o}$, we have:

The result above also holds for the Oja update up to absolute constants.

6 Related work

There is an extensive line of work analyzing PCA from the statistical perspective, in which the convergence of various estimators is characterized under certain conditions, including generative models of the data and various assumptions on the covariance matrix spectrum and eigenvalue spacing . Such works do provide finite-sample guarantees, but they apply only to the batch case and/or are computationally intensive, rather than considering an efficient incremental algorithm.

Among incremental algorithms, the work of Warmuth and Kuzmin describes and analyzes worst-case online PCA, using an experts-setting algorithm with a super-quadratic per-iteration cost. More efficient general-purpose incremental PCA algorithms have lacked finite-sample analyses . There have been recent attempts to remedy this situation by relaxing the nonconvexity inherent in the problem or making generative assumptions . The present paper directly analyzes the oldest known incremental PCA algorithms under relatively mild assumptions.

Outline of proof

We now sketch the proof of Theorem 1.1; almost all the details are relegated to the appendix.

Recall that for $n\geq n_{o}$, we take $\mathcal{F}_{n}$ to be the sigma-field of all outcomes up to and including time $n$, that is, $\mathcal{F}_{n}=\sigma(V_{n_{o}},X_{n_{o}+1},\ldots,X_{n})$.

We first bound the expected improvement in $\Psi_{n}$ in each step of the Krasulina or Oja algorithms.

For any $n>n_{o}$, we can write $\Psi_{n}\leq\Psi_{n-1}+\beta_{n}-Z_{n}$, where

and where $Z_{n}$ is a ${\mathcal{F}}_{n}$-measurable random variable with the following properties:

The theorem follows from Lemmas A.4 and A.5 in the appendix. Its characterization of the two estimators is almost identical, and for simplicity we will henceforth deal only with Krasulina’s estimator. All the subsequent results hold also for Oja’s method, up to constants.

We know from Theorem 2.1 that $\Psi_{n}\leq\Psi_{n-1}+\beta_{n}-Z_{n}$, where $\beta_{n}$ is non-stochastic and $Z_{n}$ is a quantity of positive expected value. Thus, in expectation, and modulo a small additive term, $\Psi_{n}$ decreases monotonically. However, the amount of decrease at the $n$th time step can be arbitrarily small when $\Psi_{n}$ is close to 1. Thus, we need to show that $\Psi_{n}$ is eventually bounded away from 1, i.e. there exists some $\epsilon_{o}>0$ and some time $n_{o}$ such that for any $n\geq n_{o}$, we have $\Psi_{n}\leq 1-\epsilon_{o}$.

To prove this, we start with a simple recurrence for the moment-generating function of $\Psi_{n}$.

Consider a filtration $({\mathcal{F}}_{n})$ and random variables $Y_{n},Z_{n}\in{\mathcal{F}}_{n}$ such that there are two sequences of nonnegative constants, $(\beta_{n})$ and $(\zeta_{n})$, for which:

Each $Z_{n}$ takes values in an interval of length $\zeta_{n}$.

This relation shows how to define a supermartingale based on $e^{tY_{n}}$, from which we can derive a large deviation bound on $Y_{n}$.

In order to apply this to the sequence $(\Psi_{n})$, we need to first calculate the moment-generating function of its starting value $\Psi_{n_{o}}$.

Putting these pieces together yields Theorem 2.2.

3 Intermediate epochs of improvement

We have seen that, for suitable $\epsilon$ and $n_{o}$, it is likely that $\Psi_{n}\leq 1-\epsilon/d$ for all $n\geq n_{o}$. We now define a series of epochs in which $1-\Psi_{n}$ successively doubles, until $\Psi_{n}$ finally drops below $1/2$.

To do this, we specify intermediate goals $(n_{o},\epsilon_{o}),(n_{1},\epsilon_{1}),(n_{2},\epsilon_{2}),\ldots,(n_{J},\epsilon_{J})$, where $n_{o}<n_{1}<\cdots<n_{J}$ and $\epsilon_{o}<\epsilon_{1}<\cdots<\epsilon_{J}=1/2$, with the intention that:

Of course, this can only hold with a certain probability.

Let $\Omega$ denote the sample space of all realizations $(v_{n_{o}},x_{n_{o}+1},x_{n_{o}+2},\ldots)$, and $P$ the probability distribution on these sequences. We will show that, for a certain choice of $\{(n_{j},\epsilon_{j})\}$, all $J+1$ constraints (2) can be met by excluding just a small portion of $\Omega$.

We consider a specific realization $\omega\in\Omega$ to be good if it satisfies (2). Call this set $\Omega^{\prime}$:

For technical reasons, we also need to look at realizations that are good up to time $n-1$. Specifically, for each $n$, define

Crucially, this is ${\mathcal{F}}_{n-1}$-measurable. Also note that $\Omega^{\prime}=\bigcap_{n>n_{o}}\Omega_{n}^{\prime}$.

Assume that $\gamma_{n}=c/n$, where $c=c_{o}/(2(\lambda_{1}-\lambda_{2}))$ and $c_{o}>0$. Pick any $0<\delta<1$ and select a schedule $(n_{o},\epsilon_{o}),\ldots,(n_{J},\epsilon_{J})$ that satisfies the conditions

as well as $n_{o}\geq(20c^{2}B^{2}/\epsilon_{o}^{2})\ln(4/\delta)$. Then $\mbox{\rm Pr}(\Omega^{\prime})\geq 1-\delta$.

The first step towards proving this theorem is bounding the moment-generating function of $\Psi_{n}$ in terms of that of $\Psi_{n-1}$.

Suppose $n>n_{j}$. Suppose also that $\gamma_{n}=c/n$, where $c=c_{o}/(2(\lambda_{1}-\lambda_{2}))$. Then for any $t>0$,

A repeated application of Lemmas 2.7 and 2.8 yields the following.

Suppose that conditions (3) hold. Then for $0\leq j<J$ and any $t>0$,

Now that we have bounds on the moment-generating functions of intermediate $\Psi_{n}$, we can apply martingale deviation bounds, as in Lemma 2.4, to obtain the following, from which Theorem 2.6 ensues.

Assume conditions (3) hold. Pick any $0<\delta<1$, and set $n_{o}\geq(20c^{2}B^{2}/\epsilon_{o}^{2})\ln(4/\delta)$. Then

4 The final epoch

Recall the definition of the intermediate goals $(n_{j},\epsilon_{j})$ in (2), (3). The final epoch is the period $n\geq n_{J}$, at which point $\Psi_{n}\leq 1/2$. The following consequence of Lemmas A.4 and 2.8 captures the rate at which $\Psi$ decreases during this phase.

where $\alpha_{n}=(\lambda_{1}-\lambda_{2})\gamma_{n}$ and $\beta_{n}=(B^{2}/4)\gamma_{n}^{2}$.

By solving this recurrence relation, and piecing together the various epochs, we get the overall convergence result of Theorem 1.1.

Note that Lemma 2.11 closely resembles the recurrence relation followed by the squared $L^{2}$ distance from the optimum of stochastic gradient descent (SGD) on a strongly convex function . As $\Psi_{n}\to 0$, the incremental PCA algorithms we study have convergence rates of the same form as SGD in this scenario.

Experiments

When performing PCA in practice with massive $d$ and a large/growing dataset, an incremental method like that of Krasulina or Oja remains practically viable, even as quadratic-time and -memory algorithms become increasingly impractical. Arora et al. have a more complete discussion of the empirical necessity of incremental PCA algorithms, including a version of Oja’s method which is shown to be extremely competitive in practice.

Since the efficiency benefits of these types of algorithms are well understood, we now instead focus on the effect of the learning rate on the performance of Oja’s algorithm (results for Krasulina’s are extremely similar). We use the CMU PIE faces , consisting of 11554 images of size $32\times 32$, as a prototypical example of a dataset with most of its variance captured by a few PCs, as shown in Fig. 1. We set $n_{0}=0$.

We expect from Theorem 1.1 and the discussion in the introduction that varying $c$ (the constant in the learning rate) will influence the overall rate of convergence. In particular, if $c$ is low, then halving it can be expected to halve the exponent of $n$, and the slope of the log-log convergence graph (ref. the remark after Thm. 1.1). This is exactly what occurs in practice, as illustrated in Fig. 2. The dotted line in that figure is a convergence rate of $1/n$, drawn as a guide.

Open problems

Several fundamental questions remain unanswered. First, the convergence rates of the two incremental schemes depend on the multiplier $c$ in the learning rate $\gamma_{n}$. If it is too low, convergence will be slower than $O(1/n)$. If it is too high, the constant in the rate of convergence will be large. Is there a simple and practical scheme for setting $c$?

where the second step orthonormalizes the columns, for instance by Gram-Schmidt. It would be interesting to characterize the rate of convergence of this scheme.

Finally, our analysis applies to a modified procedure in which the starting time $n_{o}$ is artificially set to a large constant. This seems unnecessary in practice, and it would be useful to extend the analysis to the case where $n_{o}=0$.

The authors are grateful to the National Science Foundation for support under grant IIS-1162581.

References

Appendix A Expected per-step change in potential

$\|\xi_{n}\|^{2}\leq B^{2}\|V_{n-1}\|^{2}/4$.

For (a), let $X_{n}^{\bot}$ denote the component of $X_{n}$ orthogonal to $V_{n-1}$. Then

For (b), note from the previous formulation that $\|\xi_{n}\|^{2}=(V_{n-1}\cdot X_{n})^{2}\|X_{n}^{\bot}\|^{2}\leq\|V_{n-1}\|^{2}\|X_{n}\|^{4}/4$.

For (d), we use $\|V_{n}\|^{2}=\|V_{n-1}+\gamma_{n}\xi_{n}\|^{2}=\|V_{n-1}\|^{2}+\gamma_{n}^{2}\|\xi_{n}\|^{2}\geq\|V_{n-1}\|^{2}$. ∎

We now check that $(V_{n}\cdot v^{*})^{2}$ grows in expectation with each iteration.

$(V_{n}\cdot v^{*})^{2}\geq(V_{n-1}\cdot v^{*})^{2}+2\gamma_{n}(V_{n-1}\cdot v^{*})(\xi_{n}\cdot v^{*})$.

Part (a) follows directly from the update rule:

In order to use Lemma A.2 to bound the change in potential $\Psi_{n}$, we need to relate $\Psi_{n}$ to the quantity $\lambda_{1}-G(V_{n})$.

For any $n\geq n_{o}$, we have $\lambda_{1}-G(V_{n})\geq(\lambda_{1}-\lambda_{2})\Psi_{n}$.

It is easiest to think of $V_{n}$ in the eigenbasis of $A$: the component of $V_{n}$ in direction $v^{*}$ is $V_{n}\cdot v^{*}$, and the orthogonal component is $V_{n}^{\bot}=V_{n}-(V_{n}\cdot v^{*})v^{*}$. Then

We can now explicitly bound the expected change in $\Psi_{n}$ in each iteration.

For any $n>n_{o}$, we can write $\Psi_{n}\leq\Psi_{n-1}+\beta_{n}-Z_{n}$, where $\beta_{n}=\gamma_{n}^{2}B^{2}/4$ and where

is a ${\mathcal{F}}_{n}$-measurable random variable with the following properties:

which is $\Psi_{n-1}+\beta_{n}-Z_{n}$. The conditional expectation of $Z_{n}$ can be determined from Lemma A.2(b):

and this can be lower-bounded using Lemma A.3.

Finally, we need to determine the range of possible values of $Z_{n}$. By expanding $\xi_{n}$, we get

Since $\|X_{n}\|^{2}\leq B$, we see that $Z_{n}$ must lie in the range $\pm 4\gamma_{n}B$. ∎

A.2 The change in potential of the Oja update

Since our bounds are on the potential function $\Psi_{n}$, which is insensitive to the length of $V_{n}$, we can skip the normalization, and instead just consider the update rule

The final bounds, as well as many of the intermediate results, are almost exactly the same as for Krasulina’s estimator. Here is the analogue of Lemma A.4.

For any $n>n_{o}$, we can write $\Psi_{n}\leq\Psi_{n-1}-Z_{n}+\beta_{n}$, where $Z_{n}$ is the same as in Lemma A.4 and $\beta_{n}=5\gamma_{n}^{2}B^{2}+2\gamma_{n}^{3}B^{3}$.

where we have used $\|X_{n}\|^{2}\leq B$. Combining these,

where the final step involves some extra algebra that we have omitted. The lemma now follows by invoking $\Psi_{n}=1-({\widehat{V}}_{n}\cdot v^{*})^{2}$. ∎

B.2 Proof of Lemma 2.4

B.3 Proof of Lemma 2.5

It is well known that $V$ can be chosen by picking $d$ values $Z=(Z_{1},\ldots,Z_{d})$ independently from the standard normal distribution and then setting $V=Z/\|Z\|$. Therefore,

where $W_{1}$ is drawn from a chi-squared distribution with $d-1$ degrees of freedom and $W_{2}$ is drawn independently from a chi-squared distribution with one degree of freedom. This characterization implies that $Y$ follows the $\mbox{Beta}((d-1)/2,1/2)$ distribution: specifically, for any $0<y<1$,

The moment-generating function of this distribution is

There isn’t a closed form for this, but an upper bound on the integral can be obtained. Assuming $d\geq 3$,

where the second step uses a change of variable $z=t(1-y)$, and the fourth uses the definition of the gamma function. To finish up, we use the inequality $\Gamma(z+1/2)\leq\sqrt{z}\,\Gamma(z)$ (Lemma B.1) to get

The following inequality is doubtless standard; we give a short proof here because we are unable to find a reference.

B.4 Proof of Theorem 2.2

To make this $\leq\epsilon/d$, it suffices to take $n_{o}\geq B^{2}c^{2}d(1+32t)/(4\epsilon)$, whereupon Lemma 2.4 yields

where the last step uses Lemma 2.5. The result follows by taking $t=d/(4\epsilon)$.

Appendix C Intermediate epochs of improvement

For any $\omega\in\Omega_{n}^{\prime}$, we have $\Psi_{n-1}(\omega)\leq 1-\epsilon_{j}$. Taking expectations over $\Omega_{n}^{\prime}$, we get the lemma.

C.2 Proof of Lemma 2.8

Let $j$ be the largest index such that $n_{j}<n$. Then

Thus the expected value of $g(\Psi_{n-1})$ over $\Omega_{n}^{\prime}$ is at most the expected value over $\Omega_{n-1}^{\prime}$.

C.3 Proof of Lemma 2.9

Define $\alpha_{n}=1-(c_{o}\epsilon_{j}/n)$ and $\xi_{n}(t)=c^{2}B^{2}t(1+32t)/4n^{2}$. By Lemmas 2.7 and 2.8, for $n>n_{j}$,

By applying these inequalities repeatedly, for $n$ shrinking to $n_{j}+1$ (and $t$ shrinking as well), we get

since $\Psi_{n_{j}}(\omega)\leq 1-\epsilon_{j}$ for all $\omega\in\Omega_{n_{j}+1}^{\prime}$. We then use the summations

To prove Lemma 2.9, we note that under conditions (3),

We have used the fact that $e^{-2x}\leq 1-x$ for $0\leq x\leq 3/4$. The rest follows by applying Lemma C.1 with $n=n_{j+1}$.

C.4 Proof of Lemma 2.10

Pick any $0<j\leq J$. We will mimic the reasoning of Theorem 2.2, being careful to define martingales only on the restricted space $\Omega_{n_{j}}^{\prime}$ and with starting time $n_{j}$. Then

To finish, we pick $t=(2/\epsilon_{o})\ln(4/\delta)$. The lower bound on $n_{o}$ is also a lower bound on $n_{j-1}$, and implies that $tc^{2}B^{2}(1+32t)/4n_{j-1}\leq t\epsilon_{o}/2$, whereupon

Appendix D The final epoch

For realizations $\omega\in\Omega_{n}^{\prime}$, we have $\Psi_{n-1}(\omega)\leq 1/2$ and thus the right-hand side of the above expression is at most $(1-\alpha_{n})\Psi_{n-1}+\beta_{n}$. Using the fact that $\Omega_{n}^{\prime}$ is ${\mathcal{F}}_{n-1}$-measurable, and taking expectations over $\Omega_{n}^{\prime}$,

as claimed. The last step uses Lemma 2.8.

D.2 Proof of Theorem 1.1

Define epochs $(n_{j},\epsilon_{j})$ that satisfy the conditions of Theorem 2.6, with $\epsilon_{J}=1/2$, and with $\epsilon_{j+1}=2\epsilon_{j}$ whenever possible. Then $J=\log_{2}1/(2\epsilon_{o})$ and

By Theorem 2.6, with probability $>1-\delta$, we have $\Psi_{n}\leq 1/2$ for all $n\geq n_{J}$. More precisely, $P(\Omega_{n}^{\prime})\geq 1-\delta$ for all $n>n_{o}$.

for $a=c_{o}/2$ and $b=c^{2}B^{2}/4$. By the $a>1$ case of Lemma D.1,

which upon further simplification yields the bound of Theorem 1.1 for $a>1$.

Consider a nonnegative sequence $(u_{t}:t\geq t_{o})$, such that for some constants $a,b>0$ and for all $t>t_{o}\geq 0$,

Then, writing the zeta function $\zeta(s)=\sum_{i=1}^{\infty}i^{-s}$,

Recursively applying the given recurrence for $u_{t}$ yields

We finish by bounding the summation of $(i+1)^{a-2}$ by a definite integral, to get: