Faster Eigenvector Computation via Shift-and-Invert Preconditioning

Introduction

On a high-level, our algorithms are based on a robust analysis of the classic idea of shift-and-invert preconditioning [Saa92], which allows us to efficiently reduce eigenvector computation to approximately solving a short sequence of well-conditioned linear systems in $\lambda\mathbf{I}-\mathbf{A}^{\top}\mathbf{A}$ for some shift parameter $\lambda\approx\lambda_{1}(\mathbf{A})$. We then apply state-of-the-art stochastic gradient methods to approximately solve these linear systems.

Typically, stochastic gradient methods are used to optimize convex functions that are given as the sum of many convex components. To solve a linear system $(\mathbf{M}^{\top}\mathbf{M})x=b$ we minimize the convex function $f(x)=\frac{1}{2}x^{\top}(\mathbf{M}^{\top}\mathbf{M})x-b^{\top}x$ with components $\psi_{i}(x)=\frac{1}{2}x^{\top}\left(m_{i}m_{i}^{\top}\right)x-\frac{1}{n}b^{\top}x$ where $m_{i}$ is the $i^{th}$ row of $\mathbf{M}$. Such an approach can be used to solve systems in $\mathbf{A}^{\top}\mathbf{A}$, however solving systems in $\mathbf{B}=\lambda\mathbf{I}-\mathbf{A}^{\top}\mathbf{A}$ requires more care. We require an analysis of SVRG that guarantees convergence even when some of our components are non-convex. We give a simple analysis for this setting, generalizing recent work in the area [SS15, CR15].

2 Our Results

Overall, our robust shifted-and-inverted power method analysis gives new understanding of this classical technique. It gives a means of obtaining provably accurate results when each iteration is implemented using fast linear system solvers with weak accuracy guarantees. In practice, this reduction between approximate linear system solving and eigenvector computation shows that optimized regression libraries can be leveraged for faster eigenvector computation in many cases. Furthermore, in theory we believe that the reduction suggests computational limits inherent in eigenvector computation as seen by the often easier-to-analyze problem of linear system solving. Indeed, in Section 7, we provide evidence that in certain regimes our statistical results are optimal.

3 Previous Work

Due to its universal applicability, eigenvector computation in the offline case is extremely well studied. Classical methods, such as the QR algorithm, take roughly $O(nd^{2})$ time to compute a full eigendecomposition. This can be accelerated to $O(nd^{\omega-1})$, where $\omega<2.373$ is the matrix multiplication constant [Wil12, LG14], however this is still prohibitively expensive for large matrices. Hence, faster iterative methods are often employed, especially when only the top eigenvector (or a few of the top eigenvectors) is desired.

The result in [Sha15c] makes an important contribution in separating input size and gap dependencies using stochastic optimization techniques. Unfortunately, the algorithm requires an approximation to the eigenvalue gap and a starting vector that has a constant dot product with the top eigenvector. In [Sha15b] the analysis is extended to a random initialization, however loses polynomial factors in $d$. Furthermore, the dependencies on the stable rank and $\epsilon$ are suboptimal – we improve them to $\operatorname{sr}(\mathbf{A})$ and $\log(1/\epsilon)$ respectively, obtaining true linear convergence.

Online Eigenvector Computation

While in the offline case the primary concern is computation time, in the online, or statistical setting, research also focuses on minimizing the number of samples that are drawn from $\mathcal{D}$ in order to achieve a given accuracy. Especially sought after are results that achieve asymptotically optimal accuracy as the sample size grows large.

A large body of work focuses on improving this simple algorithm, under a variety of assumptions on $\mathcal{D}$. A common focus is on obtaining streaming algorithms, in which the storage space is just $O(d)$ - proportional to the size of a single sample. In Table 2 we give a sampling of results in this area. All listed results rely on distributional assumptions at least as strong as those given above.

The bounds given for the simple matrix Bernstein based algorithm described above, Krasulina/Oja’s Algorithm [BDF13], and SGD [Sha15a] require no additional assumptions, aside from those given at the beginning of this section. The streaming results cited for [MCJ13] and [HP14] assume $a$ is generated from a Gaussian spike model, where $a_{i}=\sqrt{\lambda_{1}}\gamma_{i}{v_{1}}+Z_{i}$ and $\gamma_{i}\sim\mathcal{N}(0,1),Z_{i}\sim\mathcal{N}(0,I_{d})$. We note that under this model, the matrix Bernstein results improve by a $\log d$ factor and so match our results in achieving asymptotically optimal convergence rate. The results of [MCJ13] and [HP14] sacrifice this optimality in order to operate under the streaming model. Our work gives the best of both works – a streaming algorithm giving asymptotically optimal results.

4 Paper Organization

Review problem definitions and parameters for our runtime and sample bounds.

Describe the shifted-and-inverted power method and show how it can be implemented using approximate system solvers.

Show how to apply SVRG to solve systems in our shifted matrix, giving our main runtime results for offline eigenvector computation.

Show how to use an online variant of SVRG to run the shifted-and-inverted power method, giving our main sampling complexity and runtime results in the statistical setting.

Show how to efficiently estimate the shift parameters required by our algorithms.

Give a lower bound in the statistical setting, showing that our results are asymptotically optimal for a wide parameter range.

Preliminaries

2 The Statistical Problem

3 Problem Parameters

We use the following additional parameters for the offline and statistical problems respectively:

Algorithmic Framework

Here we develop our robust shift-and-invert framework. In Section 3.1 we provide a basic overview of the framework and in Section 3.2 we introduce the potential function we use to measure progress of our algorithms. In Section 3.3 we show how to analyze the framework given access to an exact linear system solver and in Section 3.4 we strengthen this analysis to work with an inexact linear system solver. Finally, in Section 3.5 we discuss initializing the framework.

2 Potential Function

Our analysis of the power method focuses on the objective of maximizing the Rayleigh quotient, $x^{\top}\mathbf{\Sigma}x$ for a unit vector $x$. Note that as the following lemma shows, this has a direct correspondence to the error in maximizing $|v_{1}^{\top}x|$:

Among all unit vectors $x$ such that $\epsilon=\lambda_{1}-x^{\top}\mathbf{\Sigma}x$, a minimizer of $\left|v_{1}^{\top}x\right|$ has the form $x=(\sqrt{1-\delta^{2}})v_{1}+\delta v_{2}$ for some $\delta$. We have

In order to track the progress of our algorithm we use a more complex potential function than just the Rayleigh quotient error, $\lambda_{1}-x^{\top}\mathbf{\Sigma}x$. Our potential function $G$ is defined for $x\neq 0$ by

where $\mathbf{P}_{v_{1}}$ and $\mathbf{P}_{v_{1}^{\perp}}$ are the projections onto $v_{1}$ and the subspace orthogonal to $v_{1}$ respectively. Equivalently, we have that:

When the Rayleigh quotient error $\epsilon=\lambda_{1}-x^{\top}\mathbf{\Sigma}x$ of $x$ is small, we can show a strong relation between $\epsilon$ and $G(x)$. We prove this in two parts. We first give a technical lemma, Lemma 2, that we will use several times for bounding the numerator of $G$. We then prove the connection in Lemma 3.

Since $\mathbf{B}=\lambda\mathbf{I}-\mathbf{\Sigma}$ and since $v_{1}$ is an eigenvector of $\mathbf{\Sigma}$ with eigenvalue $\lambda_{1}$ we have

Since $v_{1}$ is an eigenvector of $\mathbf{B}$, we can write $G(x)^{2}=\frac{x^{\top}\mathbf{B}x-(v_{1}^{\top}\mathbf{B}x)(v_{1}^{\top}x)}{(v_{1}^{\top}\mathbf{B}x)(v_{1}^{\top}x)}$. Lemmas 1 and 2 then give us:

3 Power Iteration

Here we show that the shifted-and-inverted power iteration in fact makes progress with respect to our objective function given an exact linear system solver for $\mathbf{B}$. Formally, we show that applying $\mathbf{B}^{-1}$ to a vector $x$ decreases the potential function $G(x)$ geometrically.

Let $x$ be a unit vector with $\left\langle x,v_{1}\right\rangle\neq 0$ and let $\widetilde{x}=\mathbf{B}^{-1}x$, i.e. the power method update of $\mathbf{B}^{-1}$ on $x$. Then, under our assumption on $\lambda$, we have:

Note that $\widetilde{x}$ may no longer be a unit vector. However, $G(\widetilde{x},v_{1})=G(c\widetilde{x},v_{1})$ for any scaling parameter $c$, so the theorem also holds for $\widetilde{x}$ scaled to have unit norm.

Writing $x$ in the eigenbasis, we have $x=\sum_{i}\alpha_{i}v_{i}$ and $\widetilde{x}=\sum_{i}\alpha_{i}\lambda_{i}\left(\mathbf{B}^{-1}\right)v_{i}$. Since $\left\langle x,v_{1}\right\rangle\neq 0$, $\alpha_{1}\neq 0$ and by the equivalent formulation of $G(x)$ given in (1):

The challenge in using the above theorem, and any traditional analysis of the shifted-and-inverted power method, is that we don’t actually have access to $\mathbf{B}^{-1}$. In the next section we show that the shifted-and-inverted power method is robust – we still make progress on our objective function even if we only approximate $\mathbf{B}^{-1}x$ using a fast linear system solver.

4 Approximate Power Iteration

We are now ready to prove our main result. We show that each iteration of the shifted-and-inverted power method makes constant factor expected progress on our potential function assuming we:

Start with a sufficiently good $x$ and an approximation of $\lambda_{1}$

Can apply $\mathbf{B}^{-1}$ approximately using a system solver such that the function error (i.e. distance to $\mathbf{B}^{-1}x$ in the $\mathbf{B}$ norm) is sufficiently small in expectation.

Can estimate Rayleigh quotients over $\mathbf{\Sigma}$ well enough to only accept updates that do not hurt progress on the objective function too much.

This third assumption is necessary since the second assumption is quite weak. An expected progress bound on the linear system solver allows, for example, the solver to occasionally return a solution that is entirely orthogonal to $v_{1}$, causing us to make unbounded backwards progress on our potential function. The third assumption allows us to reject possibly harmful updates and ensure that we still make progress in expectation. In the offline setting, we can access $\mathbf{A}$ and are able to compute Rayleigh quotients exactly in time $\operatorname{nnz}(\mathbf{A})$ time. However, we only assume the ability to estimate quotients since in the online setting we only have access to $\mathbf{\Sigma}$ through samples from $\mathcal{D}$.

Our general theorem for the approximate power iteration, Theorem 5, assumes that we can solve linear systems to some absolute accuracy in expectation. This is not completely standard. Typically, system solver analysis assumes an initial approximation to $\mathbf{B}^{-1}x$ and then shows a relative progress bound – that the quality of the initial approximation is improved geometrically in each iteration of the algorithm. In Corollary 6 we show how to find a coarse initial approximation to $\mathbf{B}^{-1}x$, in fact just approximating $\mathbf{B}^{-1}$ with $\frac{1}{x^{\top}\mathbf{B}x}x$. Using this approximation, we show that Theorem 5 actually implies that traditional system solver relative progress bounds suffice.

Note that in both claims we measure error of the linear system solver using $\left\|\cdot\right\|_{\mathbf{B}}$. This is a natural norm in which geometric convergence is shown for many linear system solvers and directly corresponds to the function error of minimizing $f(w)=\frac{1}{2}w^{\top}\mathbf{B}w-w^{\top}x$ to compute $\mathbf{B}^{-1}x$.

$G(\widetilde{x})\leq\frac{1}{\sqrt{10}}$ and

That is, not only do we decrease our potential function by a constant factor in expectation, but we are guaranteed that the potential function will never increase beyond $1/\sqrt{10}$.

The first claim follows directly from our choice of $\widetilde{x}$ from $x$ and $\widehat{x}$. If $\widetilde{x}=x$, it holds trivially by our assumption that $G(x)\leq\frac{1}{\sqrt{10}}$. Otherwise, $\widetilde{x}=\widehat{x}$ and we know that

and by Theorem 4 and the definition of $G$ we have

Taking expectations, using that $\left|\left\langle x,v_{1}\right\rangle\right|\leq 1$, and combining these three inequalities yields

So, conditioning on making an update and changing $x$ (i.e. $\mathcal{F}$ occurring), we see that our potential function changes exactly as in the exact case (Theorem 4) with additional additive error due to our inexact linear system solve.

which then implies by Markov inequality that

Let us now show that $\mathcal{G}\subseteq\mathcal{F}$. Suppose $\mathcal{G}$ is occurs. We can bound $\left\|\widehat{x}\right\|_{2}$ as follows:

where we use Lemmas 2 and 3 to conclude that $\left|\alpha_{1}\right|\geq\sqrt{1-\frac{1}{10}}$. We now turn to showing the Rayleigh quotient condition required by $\mathcal{F}$. In order to do this, we first bound $\widehat{x}^{\top}\mathbf{B}\widehat{x}-\left(v_{1}^{\top}\mathbf{B}\widehat{x}\right)\left(v_{1}^{\top}\widehat{x}\right)$ and then use Lemma 2. We have:

Combining (4) and (5) shows that $\mathcal{G}\subseteq\mathcal{F}$ there by proving (3).

Since $\mathbf{B}$ is PSD we see that if we let $f(w)=\frac{1}{2}w^{\top}\mathbf{B}w-w^{\top}x$, then the minimizer is $\mathbf{B}^{-1}x$. Furthermore note that $\frac{1}{x^{\top}\mathbf{B}x}=\operatorname*{arg\,min}_{\beta}f(\beta x)$ and therefore

which with Theorem 5 then completes the proof. ∎

5 Initialization

Theorem 5 and Corollary 6 show that, given a good enough approximation to $v_{1}$, we can rapidly refine this approximation by applying the shifted-and-inverted power method. In this section, we cover initialization. That is, how to obtain a good enough approximation to apply these results.

We first give a simple bound on the quality of a randomly chosen start vector $x_{0}$.

Suppose $x\sim\mathcal{N}(0,\mathbf{I})$, and we initialize $x_{0}$ as $\frac{x}{\left\|x\right\|_{2}}$, then with probability greater than $1-O\left(\frac{1}{d^{10}}\right)$, we have:

where $\kappa(\mathbf{B}^{-1})=\lambda_{1}(\mathbf{B}^{-1})/\lambda_{d}(\mathbf{B}^{-1)}$.

We now show that we can rapidly decrease our initial error to obtain the required $G(x)\leq\frac{1}{\sqrt{10}}$ bound for Theorem 5.

where $\kappa(\mathbf{B}^{-1})=\lambda_{1}(\mathbf{B}^{-1})/\lambda_{d}(\mathbf{B}^{-1)}$. Then the following procedure,

after $T=O\left(\log d+\log\kappa(\mathbf{B}^{-1}))\right)$ iterations satisfies:

with probability greater than $1-O(\frac{1}{d^{10}})$.

As before, we first bound the numerator and denominator of $G(\widehat{x})$ more carefully as follows:

We now use the above estimates to bound $G(\widehat{x})$.

By Lemma 7, we know with at least probability $1-O(\frac{1}{d^{10}})$, we have $G(x_{0})\leq\sqrt{\kappa(\mathbf{B}^{-1})}d^{10.5}$.

Conditioned on high probability result of $G(x_{0})$, we now use induction to prove $G(x_{t})\leq G(x_{0})$. It trivially holds for $t=0$. Suppose we now have $G(x)\leq G(x_{0})$, then by the condition in Theorem 8 and Markov inequality, we know with probability greater than $1-\frac{1}{100\sqrt{\kappa(\mathbf{B}^{-1})}d^{10.5}}$ we have:

The last inequality uses Corollary 6 with the fact that $\lambda_{2}\left(\mathbf{B}^{-1}\right)\leq\frac{1}{100}\lambda_{1}\left(\mathbf{B}^{-1}\right)$. Therefore, we have: We will have:

Finally, by union bound, we know with probability greater than $1-O(\frac{1}{d^{10}})$ in $T=O(\log d+\log\kappa(\mathbf{B}^{-1}))$ steps, we have:

Offline Eigenvector Computation

In this section we show how to instantiate the framework of Section 3 in order to compute an approximate top eigenvector in the offline setting. As discussed, in the offline setting we can trivially compute the Rayleigh quotient of a vector in $\operatorname{nnz}(\mathbf{A})$ time as we have explicit access to $\mathbf{A}^{\top}\mathbf{A}$. Consequently the bulk of our work in this section is to show how we can solve linear systems in $\mathbf{B}$ efficiently in expectation, allowing us to apply Corollary 6 of Theorem 5.

In Section 4.1 we first show how Stochastic Variance Reduced Gradient (SVRG) [JZ13] can be adapted to solve linear systems of the form $\mathbf{B}x=b$. If we wanted, for example, to solve a linear system in a positive definite matrix like $\mathbf{A}^{\top}\mathbf{A}$, we would optimize the objective function $f(x)=\frac{1}{2}x^{\top}\mathbf{A}^{\top}\mathbf{A}x-b^{\top}x$. This function can be written as the sum of $n$ convex components, $\psi_{i}(x)=\frac{1}{2}x^{\top}\left(a_{i}a_{i}^{\top}\right)x-\frac{1}{n}b^{\top}x$. In each iteration of traditional gradient descent, one computes the full gradient of $f(x_{i})$ and takes a step in that direction. In stochastic gradient methods, at each iteration, a single component is sampled, and the step direction is based only on the gradient of the sampled component. Hence, we avoid a full gradient computation at each iteration, leading to runtime gains.

Unfortunately, while we have access to the rows of $\mathbf{A}$ and so can solve systems in $\mathbf{A}^{\top}\mathbf{A}$, it is less clear how to solve systems in $\mathbf{B}=\lambda\mathbf{I}-\mathbf{A}^{\top}\mathbf{A}$. To do this, we will split our function into components of the form $\psi_{i}(x)=\frac{1}{2}x^{\top}\left(w_{i}\mathbf{I}-a_{i}a_{i}^{\top}\right)x-\frac{1}{n}b^{\top}x$ for some set of weights $w_{i}$ with $\sum_{i\in[n]}w_{i}=\lambda$.

Importantly, $(w_{i}\mathbf{I}-a_{i}a_{i}^{\top})$ may not be positive semidefinite. That is, we are minimizing a sum of functions which is convex, but consists of non-convex components. While recent results for minimizing such functions could be applied directly [SS15, CR15] here we show how to obtain stronger results by using a more general form of SVRG and analyzing the specific properties of our function (i.e. the variance).

With our solvers in place, in Section 4.3 we pull our results together, showing how to use these solvers in the framework of Section 3 to give faster running times for offline eigenvector computation.

Here we provide a sampling based algorithm for solving linear systems in $\mathbf{B}$. In particular we provide an algorithm for solving the more general problem where we are given a strongly convex function that is a sum of possibly non-convex functions that obey smoothness properties. We provide a general result on bounding the progress of an algorithm that solves such a problem by non-uniform sampling in Theorem 9 and then in the remainder of this section we show how to bound the requisite quantities for solving linear systems in $\mathbf{B}$.

where $\overline{S}$ is a variance parameter, then for all $m\geq 1$ we have

Consequently, if we pick $\eta$ to be a sufficiently small multiple of $1/\bar{S}$ then when $m=O(\overline{S}/\mu)$ we can decrease the error by a constant multiplicative factor in expectation.

We now apply the fact that $\left\|x+y\right\|_{2}^{2}\leq 2\left\|x\right\|_{2}^{2}+2\left\|y\right\|_{2}^{2}$ to give:

And summing over all iterations and taking expectations we have:

Theorem 9 immediately yields a solver for $\mathbf{B}x=b$. Finding the minimum norm solution to this system is equivalent to minimizing $f(x)=\frac{1}{2}x^{\top}\mathbf{B}x-b^{\top}x$. If we take the common approach of applying a smoothness bound for each $\psi_{i}$ along with a strong convexity bound on $f(x)$ we obtain:

so we have $\sum_{i\in[n]}\psi_{i}(x)=f(x)=\frac{1}{2}x^{\top}\mathbf{B}x-b^{\top}x$. Setting $p_{i}=\frac{\left\|a_{i}\right\|_{2}^{2}}{\left\|\mathbf{A}\right\|_{F}^{2}}$ for all $i$, we have

where the last step uses that $\lambda\leq 2\lambda_{1}\leq 2\left\|\mathbf{A}\right\|_{F}^{2}$ so $\frac{\lambda}{\left\|\mathbf{A}\right\|_{F}^{2}}\leq 2$. ∎

(Improved Variance Bound for SVRG) For $i\in[n]$ let

so we have $\sum_{i\in[n]}\psi_{i}(x)=f(x)=\frac{1}{2}x^{\top}\mathbf{B}x-b^{\top}x$. Setting $p_{i}=\frac{\left\|a_{i}\right\|_{2}^{2}}{\left\|\mathbf{A}\right\|_{F}^{2}}$ for all $i$, we have for all $x$

Using the gradient computation in (8) we have

Plugging the bound in Lemma 11 into Theorem 9 we have:

The procedure requires $O\left(\operatorname{nnz}(\mathbf{A})\right)$ time to initially compute $\bigtriangledown f(x_{0})$, along with each $p_{i}$ and the step size $\eta$ which depend on $\left\|\mathbf{A}\right\|_{F}^{2}$ and the row norms of $\mathbf{A}$. Each iteration then just requires $O(d)$ time to compute $\bigtriangledown\psi_{i}(\cdot)$ and perform the necessary vector operations. Since there are at most $[64\overline{S}/\mu]=O\left(\frac{\lambda_{1}\left\|\mathbf{A}\right\|_{F}^{2}}{(\lambda-\lambda_{1})^{2}}\right)$ iterations, our total runtime is

2 Accelerated Solver

Theorem 12 gives a linear solver for $\mathbf{B}$ that makes progress in expectation and which we can plug into Theorems 5 and 8. However, we first show that the runtime in Theorem 12 can be accelerated in some cases. We apply a result of [FGKS15b], which shows that, given a solver for a regularized version of a convex function $f(x)$, we can produce a fast solver for $f(x)$ itself. Specifically:

in time $\mathcal{T}_{c}$. Then given any $x_{0}$, $c>0$, $\gamma>2\mu$, we can compute $x_{1}$ such that

in time $O\left(\mathcal{T}_{4\left(\frac{2\gamma+\mu}{\mu}\right)^{3/2}}\sqrt{\lceil\gamma/\mu\rceil}\log c\right).$

We first give a new variance bound on solving systems in $\mathbf{B}$ when a regularizer is used. The proof of this bound is very close to the proof given for the unregularized problem in Lemma 11.

so we have $\sum_{i\in[n]}\psi_{i}(x)=f_{\gamma,x_{0}}(x)=\frac{1}{2}x^{\top}\mathbf{B}x-b^{\top}x+\frac{\gamma}{2}\left\|x-x_{0}\right\|_{2}^{2}$. Setting $p_{i}=\frac{\left\|a_{i}\right\|_{2}^{2}}{\left\|\mathbf{A}\right\|_{F}^{2}}$ for all $i$, we have for all $x$

For simplicity we now just use the fact that $\left\|x+y\right\|_{2}^{2}\leq 2\left\|x\right\|_{2}^{2}+2\left\|y\right\|_{2}^{2}$ and apply our bound from equation (9) to obtain:

Now, $f_{\gamma,x_{0}}(\cdot)$ is $\lambda-\lambda_{1}+\gamma$ strongly convex, so

Following Theorem 12, the variance bound of Lemma 14 means that we can make constant progress in minimizing $f_{\gamma,x_{0}}(x)$ in $O\left(\operatorname{nnz}(\mathbf{A})+dm\right)$ time where $m=O\left(\frac{\gamma^{2}+12\lambda_{1}\left\|\mathbf{\Sigma}\right\|_{F}^{2}}{(\lambda-\lambda_{1}+\gamma)^{2}}\right)$. So, for $\gamma\geq 2(\lambda-\lambda_{1})$ we can make $4\left(\frac{2\gamma+(\lambda-\lambda_{1})}{\lambda-\lambda_{1}}\right)^{3/2}$ progress, as required by Lemma 13 in time $O\left(\left(\operatorname{nnz}(\mathbf{A})+dm\right)\cdot\log\left(\frac{\gamma}{\lambda-\lambda_{1}}\right)\right)$ time. Hence by Lemma 13 we can make constant factor expected progress in minimizing $f(x)$ in time:

3 Shifted-and-Inverted Power Method

Finally, we are able to combine the solvers from Sections 4.1 and 4.2 with the framework of Section 3 to obtain faster algorithms for top eigenvector computation.

Online Eigenvector Computation

Here we show how to apply the shifted-and-inverted power method framework of Section 3 to the online setting. This setting is more difficult than the offline case. As there is no canonical matrix $\mathbf{A}$, and we only have access to the distribution $\mathcal{D}$ through samples, in order to apply Theorem 5 we must show how to both estimate the Rayleigh quotient (Section 5.1) as well as solve the requisite linear systems in expectation (Section 5.2).

After laying this ground work, our main result is given in Section 5.3. Ultimately, the results in this section allow us to achieve more efficient algorithms for computing the top eigenvector in the statistical setting as well as improve upon the previous best known sample complexity for top eigenvector computation. As we show in Section 7 the bounds we provide in this section are in fact tight for general distributions.

Here we show how to estimate the Rayleigh quotient of a vector with respect to $\mathbf{\Sigma}$. Our analysis is standard – we first approximate the Rayleigh quotient by its empirical value on a batch of $k$ samples and prove using Chebyshev’s inequality that the error on this sample is small with constant probability. We then repeat this procedure $O(\log(1/p))$ times and output the median. By a Chernoff bound this yields a good estimate with probability $1-p$. The formal statement of this result and its proof comprise the remainder of this subsection.

Given $\epsilon\in(0,1]$, $p\in$, and unit vector $x$ set $k=\lceil 4\operatorname{v}(\mathcal{D})\epsilon^{-2}\rceil$ and $m=O(\log(1/p))$. For all $i\in[k]$ and $j\in[m]$ let $a_{i}^{(j)}$ be drawn independently from $\mathcal{D}$ and set $R_{i,j}=x^{\top}a_{i}^{(j)}(a_{i}^{(j)})^{\top}x$ and $R_{j}=\frac{1}{k}\sum_{i\in[k]}R_{i,j}$. If we let $z$ be median value of the $R_{j}$ then with probability $1-p$ we have $\left|z-x^{\top}\mathbf{\Sigma}x\right|\leq\epsilon\lambda_{1}$.

The median $z$ satisfies $|z-x^{\top}\mathbf{\Sigma}x|\leq\epsilon$ as more than half of the $R_{j}$ satisfy $|R_{j}-x^{\top}\mathbf{\Sigma}x|\leq\epsilon$. This happens with probability $1-p$ by Chernoff bound, our choice of $m$ and (12). ∎

2 Solving the Linear system

The performance of streaming SVRG [FGKS15a] is governed by three regularity parameters. As in the offline case, we use the fact that $f(\cdot)$ is $\mu$-strongly convexity for $\mu=\lambda-\lambda_{1}$ and we require a smoothness parameter, denoted $\overline{S}$, that satisfies:

Furthermore, we require an upper bound the variance, denoted $\sigma^{2}$, that satisfies:

With the following two lemmas we bound these parameters.

Our proof is similar to the one for Lemma 10.

Furthermore, since $\mathbf{B}^{-1}\preceq\frac{1}{\lambda-\lambda_{1}}\mathbf{I}$ we have

Combining these three equations yields the result. ∎

With the regularity parameters bounded we can apply the streaming SVRG algorithm of [FGKS15a] to solve systems in $\mathbf{B}$. We encapsulate the core iterative step of Algorithm $1$ of [FGKS15a] as follows:

and return $x_{\widetilde{m}}$ as the output.

The accuracy of the above iterative step is proven in Theorem 4.1 of [FGKS15a], which we include, using our notation below:

Using Theorem 22 we can immediately obtain the following guarantee for solve system in $\mathbf{B}$:

Using the inequality $(x+y)^{2}\leq 2x^{2}+2y^{2}$ we have that

Now the number of samples used to compute $x$ is clearly at most $m+k$ Now

3 Online Shifted-and-Inverted Power Method

We now apply the results in Section 5.1 and Section 5.2 to the shifted-and-inverted power method framework of Section 3 to give our main result in the online setting, an algorithm that quickly refines a coarse approximation to $v_{1}$ into a finer approximation.

Parameter Estimation for Offline Eigenvector Computation

In this section, for simplicity we initially assume that we have oracle access to compute $\mathbf{B}_{\lambda}^{-1}x$ for any given $x$, and any $\lambda>\lambda_{1}$. We will then show how to achieve the same results when we can only compute $\mathbf{B}_{\lambda}^{-1}x$ approximately. We use a result of [MM15] that gives gap free bounds for computing eigenvalues using the power method. The following is a specialization of Theorem 1 from [MM15]:

Throughout the proof, we assume $\alpha$ is picked to be some large constant - e.g. $\alpha>100$. Theorem 26 implies:

Conditioning on the event that Theorem 26 holds for all iterates $i$, then the iterates of Algorithm 1 satisfy:

The proof can be decomposed into two parts:

Part I (Lines 3-4): Theorem 26 tells us that $\widetilde{\lambda}_{1}^{\left(0\right)}\geq\left(1-\frac{1}{\alpha}\right)\lambda_{1}$. This means that we have

Part II (Lines 5-6): Consider now iteration $i$. We now apply Theorem 26 to the matrix $\left(\overline{\lambda}^{\left(i-1\right)}\mathbf{I}-\mathbf{A}^{T}\mathbf{A}\right)^{-1}$. The top eigenvalue of this matrix is $\left(\overline{\lambda}^{\left(i-1\right)}-\lambda_{1}\right)^{-1}$. This means that we have $\left(1-\frac{1}{\alpha}\right)\left(\overline{\lambda}^{\left(i-1\right)}-\lambda_{1}\right)^{-1}\leq\widehat{\lambda}_{1}^{\left(i\right)}\leq\left(\overline{\lambda}^{\left(i-1\right)}-\lambda_{1}\right)^{-1}$, and hence we have,

Since $\left(\overline{\lambda}^{\left(i-1\right)}-\lambda_{2}\right)^{-1}$ is the second eigenvalue of the matrix $\left(\overline{\lambda}^{\left(i-1\right)}\mathbf{I}-\mathbf{A}^{T}\mathbf{A}\right)^{-1}$, Theorem 26 tells us that

This immediately yields the first claim. For the second claim, we notice that

where $\left(\zeta_{1}\right)$ follows from the first claim of this lemma, and $\left(\zeta_{2}\right)$ follows from Lemma 27. ∎

We now state and prove the main result in this section:

This means that the exit condition on Line $6$ must be triggered in $\overline{i}+1$ iteration, proving the first part of the lemma.

For upper bound, by Lemmas 27, 28 and exit condition we know:

Note that, although we proved the upper bound and lower bound in Theorem 29 with specific constants coefficient $\frac{1}{8}$ and $\frac{1}{120}$, this analysis can easily be extended to any smaller constants by modifying the constant in the exit condition, and choosing $\alpha$ larger. Also in the failure probability

Finally, we can also bound the runtime of algorithm 1, when we use SVRG based approximate linear system solvers for $\mathbf{B}_{\lambda}$.

Lower Bounds

Here we show that our online eigenvector estimation algorithm (Theorem 25) is asymptotically optimal - as sample size grows large it achieves optimal accuracy as a function of sample size. We rely on the following lower bound for eigenvector estimation in the Gaussian spike model:

where $\iota_{i}\sim\mathcal{N}(0,1)$, and $Z_{i}\sim\mathcal{N}(0,I_{d})$. Let $\hat{v}$ be some estimator of the top eigenvector $v^{\star}$. Then, there is some universal constant $c_{0}$, so that for $n$ sufficiently large, we have:

Suppose the claim of theorem is not true, then there exist some estimator $\hat{v}$ so that

holds for all distribution $\mathcal{D}$, and for any fixed constant $c^{\prime}$ when $n$ is sufficiently large.

Let distribution $\mathcal{D}$ be the Gaussian Spike Model specified by Eq.(16), then by calculation, it’s not hard to verify that:

Gap-Free Bounds

Let $\epsilon$ be our error parameter and $m$ be the number of eigenvalues of $\Sigma$ that are $\geq(1-\epsilon/2)\lambda_{1}$. Choose $\lambda=\lambda_{1}+\epsilon/100$. We have $\lambda_{1}(\mathbf{B}^{-1})=\frac{100}{\epsilon\lambda_{1}}$. For $i>m$ we have $\lambda_{i}(\mathbf{B}^{-1})<\frac{2}{\epsilon\lambda_{1}}$. $\kappa(\mathbf{B}^{-1})\leq\frac{100}{\epsilon}$.

Let $\mathbf{V}_{b}$ have columns equal to all bottom eigenvectors with eigenvalues $\lambda_{i}<(1-\epsilon/2)\lambda_{1}$. Let $\mathbf{V}_{t}$ have columns equal to the $m$ remaining top eigenvectors. We define a simple modified potential:

We have the following Lemma connecting this potential function to eigenvalue error:

For unit $x$, if $\bar{G}(x)\leq c\sqrt{\epsilon}$ for sufficiently small constant $c$ then $\lambda_{1}-x^{\top}\mathbf{\Sigma}x\leq\epsilon\lambda_{1}$.

So if $\bar{G}(x)\leq c\sqrt{\epsilon}$ then $\left\|\mathbf{P}_{\mathbf{V}_{t}}x\right\|^{2}_{2}c^{2}\epsilon\geq\left\|\mathbf{P}_{\mathbf{V}_{b}}x\right\|^{2}_{2}$ and since $\left\|\mathbf{P}_{\mathbf{V}_{t}}x\right\|^{2}_{2}+\left\|\mathbf{P}_{\mathbf{V}_{b}}x\right\|^{2}_{2}=1$, this gives $\left\|\mathbf{P}_{\mathbf{V}_{t}}x\right\|^{2}_{2}\geq\frac{1}{1+c^{2}\epsilon}$. So we have $x^{T}\mathbf{\Sigma}x\geq\mathbf{P}_{\mathbf{V}_{t}}x^{T}\mathbf{\Sigma}x\mathbf{P}_{\mathbf{V}_{t}}\geq\frac{(1-\epsilon/2)\lambda_{1}}{1+c^{2}\epsilon}\geq 1-\epsilon$ for small enough $c$, giving the lemma. ∎

We now follow the proof of Lemma 8, which is actually much simpler in the gap-free case.

after $T=O\left(\log d/\epsilon\right)$ iterations satisfies:

with probability greater than $1-O(\frac{1}{d^{10}})$.

By Lemma 7, we know with at least probability $1-O(\frac{1}{d^{10}})$, we have $\bar{G}(x_{0})\leq G(x_{0})\leq\sqrt{\kappa(\mathbf{B}^{-1})}d^{10.5}=\frac{100d^{10.5}}{\epsilon}$. We want to show by induction that at iteration $i$ we have $\bar{G}(x_{i})\leq\frac{1}{2^{i}}\cdot\frac{100d^{10.5}}{\epsilon}$, which will give us the lemma if we set $T=\log_{2}\left(\frac{100d^{10.5}}{c\epsilon^{1.5}}\right)=O(\log(d/\epsilon))$.

Initially, we have with high probability, by the argument in Lemma 7, $\alpha_{1}\geq\frac{1}{d^{10}}$ so we have $\left\|\mathbf{P}_{v_{1}}\left(\widehat{x}\right)\right\|_{\mathbf{B}}\geq\frac{\lambda_{1}(\mathbf{B}^{-1})}{2}\sqrt{\frac{\alpha_{1}^{2}}{\lambda_{1}(\mathbf{B}^{-1})}}$. This also holds by induction in each iteration.

Let $\hat{\alpha}_{1}=|v_{1}^{\top}\widehat{x}|/\left\|\widehat{x}\right\|_{2}$. $\left\|\mathbf{P}_{v_{1}}\left(\widehat{x}\right)\right\|_{\mathbf{B}}^{2}=\frac{\hat{\alpha}_{1}^{2}\left\|\hat{x}\right\|_{2}^{2}}{\lambda_{1}(\mathbf{B}^{-1})}$ so we have

and since $\left\|\hat{x}\right\|_{2}^{2}\leq 2\left(\left\|\mathbf{B}^{-1}x\right\|_{2}^{2}+2\left\|\xi\right\|_{2}^{2}\right)\leq\lambda_{1}(\mathbf{B}^{-1})^{2}+2\frac{\epsilon^{6}}{(3000d^{21})^{2}}\leq\lambda_{1}(\mathbf{B}^{-1})^{2}\left(2+2\frac{\epsilon^{6}}{(3000d^{21})^{2}}\right)$ we have:

So over all $\log_{2}\left(\frac{100d^{10.5}}{c\epsilon^{1.5}}\right)$ iterations, we always have $\hat{\alpha}_{1}^{2}\geq\frac{1}{d^{10}}\cdot\left(\frac{c\epsilon^{1.5}}{100d^{10.5}}\right)^{\log_{2}3}$ and so $\frac{\epsilon^{6}}{(3000d^{21})^{2}}<<1/2\alpha_{1}^{2}$. Combining the above bounds:

Finally, we combine Theorem 34 with the SVRG based solvers of Theorem 12 and 15 to obtain:

Let $\mathbf{B}=\lambda\mathbf{I}-\mathbf{A}^{\top}\mathbf{A}$ for $\lambda=\left(1+\frac{\epsilon}{100}\right)$ and let $x_{0}\sim\mathcal{N}(0,\mathbf{I})$ be a random initial vector. Running the inverted power method on $\mathbf{B}$ initialized with $x_{0}$, using the SVRG solver from Theorem 12 to approximately apply $\mathbf{B}^{-1}$ at each step, returns $x$ such that with probability $1-O\left(\frac{1}{d^{10}}\right)$, $x^{\top}\mathbf{\Sigma}x\geq(1-\epsilon)\lambda_{1}$ in time

Let $\mathbf{B}=\lambda\mathbf{I}-\mathbf{A}^{\top}\mathbf{A}$ for $\lambda=\left(1+\frac{\epsilon}{100}\right)$ and let $x_{0}\sim\mathcal{N}(0,\mathbf{I})$ be a random initial vector. Running the inverted power method on $\mathbf{B}$ initialized with $x_{0}$, using the SVRG solver from Theorem 15 to approximately apply $\mathbf{B}^{-1}$ at each step, returns $x$ such that with probability $1-O\left(\frac{1}{d^{10}}\right)$, $x^{\top}\mathbf{\Sigma}x\geq(1-\epsilon)\lambda_{1}$ in total time

Acknowledgements

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

References

Appendix A Appendix

We can any unit vector $y$ as $y=c_{1}v_{1}+c_{2}v_{2}$ where $v_{2}$ is the component of $x$ orthogonal to $v_{1}$ and $c_{1}^{2}+c_{2}^{2}=1$. We know that

We have $x^{\top}\mathbf{B}\mathbf{B}^{\top}x=c_{1}^{2}(v_{1}^{\top}\mathbf{B}^{\top}\mathbf{B}v_{1})+c_{2}^{2}(v_{2}^{\top}\mathbf{B}^{\top}\mathbf{B}v_{2})+2c_{1}c_{2}\cdot v_{2}^{\top}\mathbf{B}^{\top}\mathbf{B}v_{1}$.

We want to bound $c_{1}\geq 1-\epsilon$ so $c_{1}^{2}\geq 1-O(\epsilon)$. Since $x$ is the top eigenvector of $\mathbf{BB}^{\top}$ we have:

This means we need have $1-c_{1}^{2}\leq O(\epsilon)$ meaning $c_{1}^{2}\geq 1-O(\epsilon)$ as desired. ∎

Let $x$ be a unit vector with $\left\langle x,v_{1}\right\rangle\neq 0$ and let $\widetilde{x}=\mathbf{B}^{-1}w$, i.e. the power method update of $\mathbf{B}^{-1}$ on $x$. Then, we have both:

(17) was already shown in Lemma 4. We show (18) similarly.

Writing $x$ in the eigenbasis of $\mathbf{B}^{-1}$, we have $x=\sum_{i}\alpha_{i}v_{i}$ and $\widetilde{x}=\sum_{i}\alpha_{i}\lambda_{i}\left(\mathbf{B}^{-1}\right)v_{i}$. Since $\left\langle x,v_{1}\right\rangle\neq 0$, $\alpha_{1}\neq 0$ and we have: