Randomized Block Krylov Methods for Stronger and Faster Approximate Singular Value Decomposition

Introduction

Among countless applications, the SVD is used for optimal low-rank approximation and principal component analysis (PCA)Typically after mean centering $\mathbf{A}$’s columns or rows, depending on which principal components we want.. Specifically, for $k<r$, a partial SVD can be used to construct a rank $k$ approximation $\mathbf{A}_{k}$ such that both $\|\mathbf{A}-\mathbf{A}_{k}\|_{F}$ and $\|\mathbf{A}-\mathbf{A}_{k}\|_{2}$ are as small as possible. We simply set $\mathbf{A}_{k}=\mathbf{U}_{k}\mathbf{U}_{k}^{T}\mathbf{A}$. That is, $\mathbf{A}_{k}$ is $\mathbf{A}$ projected onto the space spanned by its top $k$ singular vectors.

For principal component analysis, $\mathbf{A}$’s top singular vector $\mathbf{u}_{1}$ provides a top principal component, which describes the direction of greatest variance within $\mathbf{A}$. The $i^{\text{th}}$ singular vector $\mathbf{u}_{i}$ provides the $i^{\text{th}}$ principal component, which is the direction of greatest variance orthogonal to all higher principal components. Formally, denoting $\mathbf{A}$’s $i^{\text{th}}$ singular value as $\sigma_{i}$,

Traditional SVD algorithms are expensive, typically running in $O(nd^{2})$ timeThis is somewhat of an oversimplicifcation. By the Abel-Ruffini Theorem, an exact SVD is incomputable even with exact arithmetic . Accordingly, all SVD algorithm are inherently iteratively. Nevertheless, traditional methods including the ubiquitous QR algorithm obtain superlinear convergence rates for the low-rank approximation problem. In any reasonable computing environment, they can be taken to run in $O(nd^{2})$ time.. Hence, there has been substantial research on randomized techniques that seek nearly optimal low-rank approximation and PCA . These methods are quickly becoming standard tools in practice and implementations are widely available , including in popular learning libraries like scikit-learn .

Recent work focuses on algorithms whose runtimes do not depend on properties of $\mathbf{A}$. In contrast, classical literature typically gives runtime bounds that depend on the gaps between $\mathbf{A}$’s singular values and become useless when these gaps are small (which is often the case in practice – see Section 8). This limitation is due to a focus on how quickly approximate singular vectors converge to the actual singular vectors of $\mathbf{A}$. When two singular vectors have nearly identical values they are difficult to distinguish, so convergence inherently depends on singular value gaps.

Only recently has a shift in approximation goal, along with an improved understanding of randomization, allowed for algorithms that avoid gap dependence and thus run provably fast for any matrix. For low-rank approximation and PCA, we only need to find a subspace that captures nearly as much variance as $\mathbf{A}$’s top singular vectors – distinguishing between two close singular values is overkill.

The fastest randomized SVD algorithms run in $O(\operatorname{nnz}(\mathbf{A}))$ timeHere $\operatorname{nnz}(\mathbf{A})$ is the number of non-zero entries in $\mathbf{A}$ and this runtime hides lower order terms., are based on non-iterative sketching methods, and return a rank $k$ matrix $\mathbf{Z}$ with orthonormal columns $\mathbf{z}_{1},\ldots,\mathbf{z}_{k}$ satisfying

Unfortunately, as emphasized in prior work , Frobenius norm error is often hopelessly insufficient, especially for data analysis and learning applications. When $\mathbf{A}$ has a “heavy-tail” of singular values, which is common for noisy data, $\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}=\sum_{i>k}\sigma_{i}^{2}$ can be huge, potentially much larger than $\mathbf{A}$’s top singular value. This renders (1) meaningless since $\mathbf{Z}$ does not need to align with any large singular vectors to obtain good multiplicative error.

To address this shortcoming, a number of papers suggest targeting spectral norm low-rank approximation error,

which is intuitively stronger. When looking for a rank $k$ approximation, $\mathbf{A}$’s top $k$ singular vectors are often considered data and the remaining tail is considered noise. A spectral norm guarantee roughly ensures that $\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}$ recovers $\mathbf{A}$ up to this noise threshold.

Our Results

Even though the algorithm has been discussed and tested for potential improvement over Simultaneous Iteration , theoretical bounds for Krylov subspace and Lanczos methods are much more limited. As highlighted in ,

“Despite decades of research on Lanczos methods, the theory for [randomized power iteration] is more complete and provides strong guarantees of excellent accuracy, whether or not there exist any gaps between the singular values.”

Our work addresses this issue, giving the first gap independent bound for a Krylov subspace method.

2 Stronger Guarantees

In addition to runtime improvements, we target a much stronger notion of approximate SVD that is needed for many applications, but for which no gap-independent analysis was known.

Specifically, as noted in , while intuitively stronger than Frobenius norm error, $(1+\epsilon)$ spectral norm low-rank approximation error does not guarantee any accuracy in $\mathbf{Z}$ for many matricesIn fact, it does not even imply $(1+\epsilon)$ Frobenius norm error.. Consider $\mathbf{A}$ with its top $k+1$ squared singular values all equal to $10$ followed by a tail of smaller singular values (e.g. $1000k$ at $1$). $\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{2}=10$ but in fact $\|\mathbf{A}-\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}\|^{2}_{2}=10$ for any rank $k$ $\mathbf{Z}$, leaving the spectral norm bound useless. At the same time, $\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}$ is large, so Frobenius error is meaningless as well. For example, any $\mathbf{Z}$ obtains $\|\mathbf{A}-\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}\|_{F}^{2}\leq(1.01)\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}$.

With this scenario in mind, it is unsurprising that low-rank approximation guarantees fail as an accuracy measure in practice. We ran a standard sketch-and-solve approximate SVD algorithm (see Section 3.1) on SNAP/amazon0302, an Amazon product co-purchasing dataset , and achieved very good low-rank approximation error in both norms for $k=30$:

However, the approximate principal components given by $\mathbf{Z}$ are of significantly lower quality than $\mathbf{A}$’s true singular vectors (see Figure 1). We saw a similar phenomenon for the popular 20 Newsgroups dataset and several others. Additionally, the potential failure of low rank approximation measures was recently raised in .

We address this issue by introducing a per vector guarantee that requires each approximate singular vector $\mathbf{z}_{1},\ldots,\mathbf{z}_{k}$ to capture nearly as much variance as the corresponding true singular vector:

The error bound (3) is very strong in that it depends on $\epsilon\sigma^{2}_{k+1}$, meaning that it is better then relative error, i.e. $\left|\mathbf{u}_{i}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{u}_{i}-\mathbf{z}_{i}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{z}_{i}\right|\leq\epsilon\sigma_{i}^{2}$, for $\mathbf{A}$’s large singular vectors. While it is reminiscent of the bounds sought in classical numerical analysis , we stress that it does not require each $\mathbf{z}_{i}$ to converge to $\mathbf{u}_{i}$ in the presence of small singular value gaps. In fact, we show that both randomized Block Krylov Iteration and our slightly modified Simultaneous Iteration algorithmFor guarantee (3) it is important that Algorithm 1 includes post-processing steps 4 and 5 rather than just returning a basis for $\mathbf{K}$, which is sufficient for the low-rank approximation guarantees. achieve (3) in gap-independent runtimes.

3 Main Result

Our contributions are summarized in Theorem 1, whose proof appears in parts as Theorems 6 and 7 in Section 5 (runtime) and Theorems 10, 11, and 12 in Section 6 (accuracy).

With high probability, Algorithms 1 and 2 find approximate singular vectors $\mathbf{Z}=\left[\mathbf{z}_{1},\ldots,\mathbf{z}_{k}\right]$ satisfying guarantees (1) and (2) for low-rank approximation and (3) for PCA. For error $\epsilon$, Algorithm 1 requires $q=O(\log d/\epsilon)$ iterations while Algorithm 2 requires $q=O(\log d/\sqrt{\epsilon})$ iterations. Excluding lower order terms, both algorithms run in time $O(\operatorname{nnz}(\mathbf{A})kq)$.

We note that, while Simultaneous Iteration was known to achieve (2) , surprisingly we are first to prove that it gives (1), a qualitatively weaker goal.

In Section 7 we use our results to give an alternative analysis of both algorithms that does depend on singular value gaps and can offer significantly faster convergence when $\mathbf{A}$ has decaying singular values. It is possible to take further advantage of this result by running Algorithms 1 and 2 with a $\mathbf{\Pi}$ that has $>k$ columns, a simple modification for accelerating either method.

Finally, Section 8 contains a number of experiments on large data problems. We justify the importance of gap independent bounds for predicting algorithm convergence and we show that Block Krylov Iteration in fact significantly outperforms the more popular Simultaneous Iteration.

4 Comparison to Classical Bounds

Decades of work has produced a variety of gap dependent bounds for power iteration and Krylov subspace methods. We refer the reader to Saad’s standard reference . Most relevant to our work are bounds for block Krylov methods with block size equal to $k$ . Roughly speaking, with randomized initialization, these results offer guarantees equivalent to our strong equation (3) for the top $k$ singular directions after:

This bound is recovered by our Section 7 results and, when the target accuracy $\epsilon$ is smaller than the relative singular value gap $(\sigma_{k}/\sigma_{k+1}-1)$, it is tighter than our gap independent results. However, as discussed in Section 8, for high dimensional data problems where $\epsilon$ is set far above machine precision, gap independent bounds more accurately predict required iteration count.

Less comparable to our results are attempts to analyze algorithms with block size smaller than $k$ . While “small block” or single vector algorithms offer runtime advantages, it is well understood that with $b$ duplicate singular values, it is impossible to recover the top $k$ singular directions with a block of size $<b$ . More generally, large singular value clusters slow convergence, so any small block algorithm must have runtime dependence on the gaps between each adjacent pair of top $k$ singular values . We believe that obtaining simpler theoretical bounds for small block methods is an interesting direction for future work.

Background and Intuition

We will start by 1) providing background on algorithms for approximate singular value decomposition and 2) giving intuition for Simultaneous Power Iteration and Block Krylov methods and justifying why they can give strong gap-independent error guarantees.

Progress on algorithms for Frobenius norm error low-rank approximation (1) has been considerable. Work in this direction dates back to the strong rank-revealing QR factorizations of Gu and Eisenstat . They give deterministic algorithms that run in approximately $O(ndk)$ time, vs. $O(nd^{2})$ for a full SVD, but only guarantee polynomial factor Frobenius norm error.

The sketch-and-solve method is very efficient – the computation of $\mathbf{A\Pi}$ is easily parallelized and, regardless, pass-efficient in a single processor setting. Furthermore, once a small compression of $\mathbf{A}$ is obtained, it can be manipulated in fast memory to find $\mathbf{Z}$. This is not typically true of $\mathbf{A}$ itself, making it difficult to directly process the original matrix at all.

2 Spectral Norm Error via Simultaneous Iteration

Unfortunately, as discussed, Frobenius norm error is often insufficient when $\mathbf{A}$ has a heavy singular value tail. Moreover, it seems an inherent limitation of sketch-and-solve methods. The noise from $\mathbf{A}$’s lower $r-k$ singular values corrupts $\mathbf{A\Pi}$, making it impossible to extract a good partial SVD if the sum of these singular values (equal to $\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}$) is too large. In other words, any error inherently depends on the size of this tail.

In order to achieve spectral norm error (2), Simultaneous Iteration must reduce this noise down to the scale of $\sigma_{k+1}=\|\mathbf{A}-\mathbf{A}_{k}\|_{2}$. It does this by working with the powered matrix $\mathbf{A}^{q}$ .For nonsymmetric matrices we work with ($\mathbf{A}\mathbf{A}^{T})^{q}\mathbf{A}$, but present the symmetric case here for simplicity. By the spectral theorem, $\mathbf{A}^{q}$ has exactly the same singular vectors as $\mathbf{A}$, but its singular values are equal to the singular values of $\mathbf{A}$ raised to the $q^{\text{th}}$ power. Powering spreads the values apart and accordingly, $\mathbf{A}^{q}$’s lower singular values are relatively much smaller than its top singular values (see Figure 2(a) for an example).

Computing $\mathbf{A}^{q}$ directly is costly, so $\mathbf{A}^{q}\mathbf{\Pi}$ is computed iteratively. We start with a random $\mathbf{\Pi}$ and repeatedly multiply by $\mathbf{A}$ on the left. Since even a rough Frobenius norm approximation for $\mathbf{A}^{q}$ suffices, $\mathbf{\Pi}$ is often chosen to have just $k$ columns. Each iteration thus takes $O(\operatorname{nnz}(\mathbf{A})k)$ time. After $\mathbf{A}^{q}\mathbf{\Pi}$ is computed, $\mathbf{Z}$ can simply be set to a basis for its column span.

To the best of our knowledge, this approach to analyzing Simultaneous Iteration without dependence on singular value gaps began with . The technique was popularized in and its analysis improved in and . gives the first bound that directly achieves (2) with $O(\log d/\epsilon)$ power iterations. All of these papers rely on an improved understanding of the benefits of starting with a randomized $\mathbf{\Pi}$, which has developed from work on the sketch-and-solve paradigm.

3 Beating Simultaneous Iteration with Krylov Methods

As mentioned, numerous papers hint at the possibility of beating Simultaneous Iteration with block Krylov methods . In particular, , and suggest and experimentally confirm the potential of a randomized variant of the Block Lanczos algorithm, which we refer to as Block Krylov Iteration (Algorithm 2). However, none of these papers give theoretical bounds on the algorithm’s performance.

The intuition behind Block Krylov Iteration matches that of many accelerated iterative methods. Simply put, there are better polynomials than $\mathbf{A}^{q}$ for denoising tail singular values. In particular, we can use a lower degree polynomial, allowing us to compute fewer powers of $\mathbf{A}$ and thus leading to an algorithm with fewer iterations. For example, an appropriately shifted $q=O(\frac{\log d}{\sqrt{\epsilon}})$ degree Chebyshev polynomial can push the tail of $\mathbf{A}$ nearly as close to zero as $\mathbf{A}^{O(\log d/\epsilon)}$, even if the long run growth of the polynomial is much lower (see Figure 2(b)).

Block Krylov Iteration takes advantage of such polynomials by working with the Krylov subspace,

from which we can construct $p_{q}(\mathbf{A})\mathbf{\Pi}$ for any polynomial $p_{q}(\cdot)$ of degree $q$.Algorithm 2 in fact only constructs odd powered terms in $\mathbf{K}$, which is sufficient for our choice of $p_{q}(x)$. Since an effective polynomial for denoising $\mathbf{A}$ must be scaled and shifted based on the value of $\sigma_{k+1}$, we cannot easily compute it directly. Instead, we argue that the very best $k$ rank approximation to $\mathbf{A}$ lying in the span of $\mathbf{K}$ at least matches the approximation achieved by projecting onto the span of $p_{q}(\mathbf{A})\mathbf{\Pi}$. Finding this best approximation will therefore give a nearly optimal low-rank approximation to $\mathbf{A}$.

Unfortunately, there’s a catch. Perhaps surprisingly, it is not clear how to efficiently compute the best spectral norm error low-rank approximation to $\mathbf{A}$ lying in a specific subspace (e.g. $\mathbf{K}$’s span) . This challenge precludes an analysis of Krylov methods parallel to the recent work on Simultaneous Iteration. Nevertheless, we show that computing the best Frobenius error low-rank approximation in the span of $\mathbf{K}$, exactly the post-processing step taken by classic Block Lanczos and our method, will give a good enough spectral norm approximation for achieving $(1+\epsilon)$ error.

4 Stronger Per Vector Error Guarantees

Achieving the per vector guarantee of (3) requires a more nuanced understanding of how Simultaneous Iteration and Block Krylov Iteration denoise the spectrum of $\mathbf{A}$. The analysis for spectral norm low-rank approximation relies on the fact that $\mathbf{A}^{q}$ (or $p_{q}(\mathbf{A})$ for Block Krylov Iteration) blows up any singular value $\geq(1+\epsilon)\sigma_{k+1}$ to much larger than any singular value $\leq\sigma_{k+1}$. This ensures that the $\mathbf{Z}$ outputted by both algorithms aligns very well with the singular vectors corresponding to these large singular values.

If $\sigma_{k}\geq(1+\epsilon)\sigma_{k+1}$, then $\mathbf{Z}$ aligns well with all top $k$ singular vectors of $\mathbf{A}$ and we get good Frobenius norm error and the per vector guarantee (3). Unfortunately, when there is a small gap between $\sigma_{k}$ and $\sigma_{k+1}$, $\mathbf{Z}$ could miss intermediate singular vectors whose values lie between $\sigma_{k+1}$ and $(1+\epsilon)\sigma_{k+1}$. This is the case where gap dependent guarantees of classical analysis break down.

However, $\mathbf{A}^{q}$ or, for Block Krylov Iteration, some $q$-degree polynomial in our Krylov subspace, also significantly separates singular values $>\sigma_{k+1}$ from those $<(1-\epsilon)\sigma_{k+1}$. Thus, each column of $\mathbf{Z}$ at least aligns with $\mathbf{A}$ nearly as well as $\mathbf{u}_{k+1}$. So, even if we miss singular values between $\sigma_{k+1}$ and $(1+\epsilon)\sigma_{k+1}$, they will be replaced with approximate singular values $>(1-\epsilon)\sigma_{k+1}$, enough for (3).

For Frobenius norm low-rank approximation, we prove that the degree to which $\mathbf{Z}$ falls outside of the span of $\mathbf{A}$’s top $k$ singular vectors depends on the number of singular values between $\sigma_{k+1}$ and $(1-\epsilon)\sigma_{k+1}$. These are the values that could be ‘swapped in’ for the true top $k$ singular values. Since their weight counts towards $\mathbf{A}$’s tail, our total loss compared to optimal is at worst $\epsilon\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$.

Preliminaries

Before proceeding to the full technical analysis, we overview required results from linear algebra, polynomial approximation, and randomized low-rank approximation.

Let $\mathbf{\Sigma}_{k}$ be $\mathbf{\Sigma}$ with all but its largest $k$ singular values zeroed out. Let $\mathbf{U}_{k}$ and $\mathbf{V}_{k}$ be $\mathbf{U}$ and $\mathbf{V}$ with all but their first $k$ columns zeroed out. For any $k$, $\mathbf{A}_{k}=\mathbf{U}\mathbf{\Sigma}_{k}\mathbf{V^{T}}=\mathbf{U}_{k}\mathbf{\Sigma}_{k}\mathbf{V}_{k}^{T}$ is the closest rank $k$ approximation to $\mathbf{A}$ for any unitarily invariant norm, including the Frobenius norm and spectral norm . The squared Frobenius norm is given by $\|\mathbf{A}\|^{2}_{F}=\sum_{i,j}\mathbf{A}_{i,j}^{2}=\operatorname{tr}(\mathbf{AA^{T}})=\sum_{i}\sigma_{i}^{2}$. The spectral norm is given by $\|\mathbf{A}\|_{2}=\sigma_{1}$.

We often work with the remainder matrix $\mathbf{A}-\mathbf{A}_{k}$ and label it $\mathbf{A}_{r\setminus k}$. Its singular value decomposition is given by $\mathbf{A}_{r\setminus k}=\mathbf{U}_{r\setminus k}\mathbf{\Sigma}_{r\setminus k}\mathbf{V}_{r\setminus k}^{T}$ where $\mathbf{U}_{r\setminus k}$, $\mathbf{\Sigma}_{r\setminus k}$, and $\mathbf{V}_{r\setminus k}^{T}$ have their first $k$ columns zeroed.

While the SVD gives a globally optimal rank $k$ approximation for $\mathbf{A}$, both Simultaneous Iteration and Block Krylov Iteration return the best $k$ rank approximation falling within some fixed subspace spanned by a basis $\mathbf{Q}$ (with rank $\geq k$). For the Frobenius norm, this simply requires projecting $\mathbf{A}$ to $\mathbf{Q}$ and taking the best rank $k$ approximation of the resulting matrix using an SVD.

This low-rank approximation can be obtained using an SVD (equivalently, eigendecomposition) of the $m\times m$ matrix $\mathbf{M}=\mathbf{Q}^{T}(\mathbf{AA}^{T})\mathbf{Q}$. Specifically, letting $\mathbf{M}=\mathbf{\bar{U}}\mathbf{\bar{\Sigma}}^{2}\mathbf{\bar{U}}^{T}$, then:

If the SVD of $\mathbf{Q}^{T}\mathbf{A}$ is given by $\mathbf{Q}^{T}\mathbf{A}=\mathbf{\bar{U}}\mathbf{\bar{\Sigma}}\mathbf{\bar{V}}^{T}$ then $\mathbf{M}=\mathbf{Q}^{T}(\mathbf{AA}^{T})\mathbf{Q}=\mathbf{\bar{U}}\mathbf{\bar{\Sigma}}^{2}\mathbf{\bar{U}}^{T}$. So $\mathbf{Q}\left(\mathbf{Q}^{T}\mathbf{A}\right)_{k}=\mathbf{Q}\mathbf{\bar{U}}_{k}\mathbf{\bar{\Sigma}}_{k}\mathbf{\bar{V}}_{k}^{T}=\mathbf{Q}\left(\mathbf{\bar{U}}_{k}\mathbf{\bar{U}}_{k}^{T}\right)\mathbf{\bar{U}}\mathbf{\bar{\Sigma}}\mathbf{\bar{V}}^{T}=\mathbf{Q}\mathbf{\bar{U}}_{k}\mathbf{\bar{U}}_{k}^{T}\mathbf{Q}^{T}\mathbf{A}$, giving the lower matrix equality. Note that $\mathbf{Q}\mathbf{\bar{U}}_{k}$ has orthonormal columns since $\mathbf{\bar{U}}_{k}^{T}\mathbf{Q}^{T}\mathbf{Q}\mathbf{\bar{U}}_{k}=\mathbf{\bar{U}}_{k}^{T}\mathbf{I}\mathbf{\bar{U}}_{k}=\mathbf{I}_{k}$.

In general, this rank $k$ approximation does not give the best spectral norm approximation to $\mathbf{A}$ falling within $\mathbf{Q}$ . A closed form solution can be obtained using the results of , which are related to Parrott’s theorem, but we do not know how to compute this solution without essentially performing an SVD of $\mathbf{A}$. It is at least simple to show that the optimal spectral norm approximation for $\mathbf{A}$ spanned by a rank $k$ basis is obtained by projecting $\mathbf{A}$ to the basis:

2 Other Linear Algebra Tools

Throughout this paper we use $span(\mathbf{M})$ to denote the column span of the matrix $\mathbf{M}$. We say that a matrix $\mathbf{Q}$ is an orthonormal basis for the column span of $\mathbf{M}$ if $\mathbf{Q}$ has orthonormal columns and $\mathbf{Q}\mathbf{Q}^{T}\mathbf{M}=\mathbf{M}$. That is, projecting the columns of $\mathbf{M}$ to $\mathbf{Q}$ fully recovers those columns. $\mathbf{QQ}^{T}$ is the orthogonal projection matrix onto the span of $\mathbf{Q}$. $(\mathbf{QQ}^{T})(\mathbf{QQ}^{T})=\mathbf{QIQ}^{T}=\mathbf{QQ}^{T}$.

If $\mathbf{M}$ and $\mathbf{N}$ have the same dimension and $\mathbf{MN^{T}}=\mathbf{0}$ then $\|\mathbf{M}+\mathbf{N}\|^{2}_{F}=\|\mathbf{M}\|^{2}_{F}+\|\mathbf{N}\|^{2}_{F}$. This matrix Pythagorean theorem follows from writing $\|\mathbf{M+N}\|^{2}_{F}=\operatorname{tr}(\mathbf{(M+N)(M+N)^{T}})$. As an example, for any orthogonal projection $\mathbf{QQ}^{T}\mathbf{A}$, $\mathbf{A}^{T}(\mathbf{I}-\mathbf{QQ}^{T})\mathbf{QQ}^{T}\mathbf{A}=\mathbf{0}$, so $\|\mathbf{A}-\mathbf{QQ}^{T}\mathbf{A}\|_{F}^{2}=\|\mathbf{A}\|_{F}^{2}-\|\mathbf{QQ}^{T}\mathbf{A}\|_{F}^{2}.$ This implies that, since $\mathbf{A}_{k}=\mathbf{U}_{k}\mathbf{U}_{k}^{T}\mathbf{A}$ minimizes $\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}$ over all rank $k$ matrices, $\mathbf{QQ}^{T}=\mathbf{U}_{k}\mathbf{U}_{k}$ maximizes $\|\mathbf{QQ}^{T}\mathbf{A}\|_{F}^{2}$ over all rank $k$ orthogonal projections.

3 Randomized Low-Rank Approximation

Our proofs build on well known sketch-based algorithms for low-rank approximation with Frobenius norm error. A short proof of the following Lemma is in Appendix A:

For analyzing block methods, results like Lemma 4 can effectively serve as a replacement for earlier random initialization analysis that applies to single vector power and Krylov methods .

4 Chebyshev Polynomials

As outlined in Section 3.3, our proof also requires polynomials to more effectively denoise the tail of $\mathbf{A}$. As is standard for Krylov subspace methods, we use a variation on the Chebyshev polynomials. The proof of the following Lemma is relegated to Appendix A.

Given a specified value $\alpha>0$, gap $\gamma\in(0,1]$, and $q\geq 1$, there exists a degree $q$ polynomial $p(x)$ such that:

$p(\left(1+\gamma)\alpha\right)=(1+\gamma)\alpha$

$p(x)\geq x$ for all $x\geq(1+\gamma)\alpha$

$|p(x)|\leq\frac{\alpha}{2^{q\sqrt{\gamma}-1}}$ for all $x\in[0,\alpha]$

Furthermore, when $q$ is odd, the polynomial only contains odd powered monomials.

Implementation and Runtimes

We first briefly discuss runtime and implementation considerations for Algorithms 1 and 2, our randomized implementations of Simultaneous Power Iteration and Block Krylov Iteration.

Algorithm 1 can be modified in a number of ways. $\mathbf{\Pi}$ can be replaced by a random sign matrix, or any matrix achieving the guarantee of Lemma 4. $\mathbf{\Pi}$ may also be chosen with $p>k$ columns. We will discuss in detail how this approach can give improved accuracy in Section 7.

In our implementation we set $\mathbf{Z}=\mathbf{Q\mathbf{\bar{U}}_{k}}$. This ensures that, for all $l\leq k$, $\mathbf{Z}_{l}$ gives the best rank $l$ Frobenius norm approximation to $\mathbf{A}$ within the span of $\mathbf{K}$ (See Lemma 2). This is necessary for achieving per vector guarantees for approximate PCA. However, if we are only interested in computing a near optimal low-rank approximation, we can simply set $\mathbf{Z}=\mathbf{Q}$. Projecting $\mathbf{A}$ to $\mathbf{Q\mathbf{\bar{U}}_{k}}$ is equivalent to projecting to $\mathbf{Q}$ as these two matrices have the same column spans.

Additionally, since powering $\mathbf{A}$ spreads its singular values, $\mathbf{K}=(\mathbf{AA}^{T})^{q}\mathbf{A\Pi}$ could be poorly conditioned. As suggested in , to improve stability we can orthonormalize $\mathbf{K}$ after every iteration (or every few iterations). This does not change $\mathbf{K}$’s column span, so it gives an equivalent algorithm in exact arithmetic, but improves conditioning significantly.

Computing $\mathbf{K}$ requires first multiplying $\mathbf{A}$ by $\mathbf{\Pi}$, which takes $O(\operatorname{nnz}(\mathbf{A})k)$ time. Computing $\left(\mathbf{A}\mathbf{A}^{T}\right)^{i}\mathbf{A\Pi}$ given $\left(\mathbf{A}\mathbf{A}^{T}\right)^{i-1}\mathbf{A\Pi}$ then takes $O(\operatorname{nnz}(\mathbf{A})k)$ time to first multiply our $(n\times k)$ matrix by $\mathbf{A}^{T}$ and then by $\mathbf{A}$. Reorthogonalizing after each iteration takes $O(nk^{2})$ time via Gram-Schmidt or Householder reflections. This gives a total runtime of $O(\operatorname{nnz}(\mathbf{A})kq+nk^{2}q)$ for computing $\mathbf{K}$.

Finding $\mathbf{Q}$ takes $O(nk^{2})$ time. Computing $\mathbf{M}$ by multiplying from left to right requires $O(nnz(\mathbf{A})k+nk^{2})$ time. $\mathbf{M}$’s SVD then requires $O(k^{3})$ time using classical techniques. Finally, multiplying $\mathbf{\bar{U}}_{k}$ by $\mathbf{Q}$ takes time $O(nk^{2})$. Setting $q=\Theta(\log d/\epsilon)$ gives the claimed runtime. ∎

2 Block Krylov Iteration

As with Simultaneous Iteration, we can replace $\mathbf{\Pi}$ with any matrix achieving the guarantee of Lemma 4 and can use $p>k$ columns to improve accuracy. $\mathbf{Q}$ can also be computed in a number of ways. In the traditional Block Lanczos algorithm, one starts by computing an orthonormal basis for $\mathbf{A\Pi}$, the first block in the Krylov subspace. Bases for subsequent blocks are computed from previous blocks using a three term recurrence that ensures $\mathbf{Q}^{T}\mathbf{AA}^{T}\mathbf{Q}$ is block tridiagonal, with $k\times k$ sized blocks . This technique can be useful if $qk$ is large, since it is faster to compute the top singular vectors of a block tridiagonal matrix. However, computing $\mathbf{Q}$ using a recurrence can introduce a number of stability issues, and additional steps may be required to ensure that the matrix remains orthogonal .

An alternative is to compute $\mathbf{K}$ explicitly and then compute $\mathbf{Q}$ using a QR decomposition. This method is used in and . It does not guarantee that $\mathbf{Q}^{T}\mathbf{AA}^{T}\mathbf{Q}$ is block tridiagonal, but helps avoid a number of stability issues. Furthermore, if $qk$ is small, taking the SVD of $\mathbf{Q}^{T}\mathbf{AA}^{T}\mathbf{Q}$ will still be fast and typically dominated by the cost of computing $\mathbf{K}$.

As with Simultaneous Iteration, we can also orthonormalize each block of $\mathbf{K}$ after it is computed, avoiding poorly conditioned blocks and giving an equivalent algorithm in exact arithmetic.

Computing $\mathbf{K}$, including block reorthogonalization, requires $O(\operatorname{nnz}(\mathbf{A})kq+nk^{2}q)$ time. The remaining steps are analogous to those in Simultaneous Iteration except somewhat more costly as we work an $k\cdot q$ dimensional rather than $k$ dimensional subspace. Finding $\mathbf{Q}$ takes $O(n(kq)^{2})$ time. Computing $\mathbf{M}$ take $O(nnz(\mathbf{A})(kq)+n(kq)^{2})$ time and its SVD then requires $O((kq)^{3})$ time. Finally, multiplying $\mathbf{\bar{U}}_{k}$ by $\mathbf{Q}$ takes time $O(nk(kq))$. Setting $q=\Theta(\log d/\sqrt{\epsilon})$ gives the claimed runtime. ∎

Error Bounds

We next prove that both Algorithms 1 and 2 return a basis $\mathbf{Z}$ that gives relative error Frobenius (1) and spectral norm (2) low-rank approximation error as well as the per vector guarantees (3).

We start with a general approximation lemma, which gives three guarantees formalizing the intuition given in Section 3. All other proofs follow nearly immediately from this lemma.

For simplicity we assume that $k\leq r=\operatorname{rank}(\mathbf{A})\leq n,d$. However, if $k>r$ it can be seen that both algorithms still return a basis satisfying the proven guarantees. We start with a definition:

Recall that $\mathbf{A}_{l}$ is the best rank $l$ approximation to $\mathbf{A}$. This error function measures how well $\mathbf{Z}_{l}\mathbf{Z}_{l}^{T}\mathbf{A}$ approximates $\mathbf{A}$ in comparison to the optimal.

Let $m$ be the number of singular values $\sigma_{i}$ of $\mathbf{A}$ with $\sigma_{i}\geq(1+\epsilon/2)\sigma_{k+1}$. Let $w$ be the number of singular values with $\frac{1}{1+\epsilon/2}\sigma_{k}\leq\sigma_{i}<\sigma_{k}$. With probability $99/100$ Algorithms 1 and 2 return $\mathbf{Z}$ satisfying:

$\forall l\leq m$, $\mathcal{E}(\mathbf{Z}_{l},\mathbf{A})\leq(\epsilon/2)\cdot\sigma_{k+1}^{2},$

$\forall l\leq k$, $\mathcal{E}(\mathbf{Z}_{l},\mathbf{A})\leq\mathcal{E}(\mathbf{Z}_{l-1},\mathbf{A})+3\epsilon\cdot\sigma_{k+1}^{2},$

$\forall l\leq k$, $\mathcal{E}(\mathbf{Z}_{l},\mathbf{A})\leq(w+1)\cdot 3\epsilon\cdot\sigma_{k+1}^{2}.$

Property 1 captures the intuition given in Section 3.2. Both algorithms return $\mathbf{Z}$ with $\mathbf{Z}_{l}$ equal to the best Frobenius norm low-rank approximation in $span(\mathbf{K})$. Since $\sigma_{1}\geq\ldots\geq\sigma_{m}\geq(1+\epsilon/2)\sigma_{k+1}$ and our polynomials separate any values above this threshold from anything below $\sigma_{k+1}$, $\mathbf{Z}$ must align very well with $\mathbf{A}$’s top $m$ singular vectors. Thus $\mathcal{E}(\mathbf{Z}_{l},\mathbf{A})$ is very small for all $l\leq m$.

Property 2 captures the intuition of Section 3.4 – outside of the largest $m$ singular values, $\mathbf{Z}$ still performs well. We may fail to distinguish between vectors with values between $\frac{1}{1+\epsilon/2}\sigma_{k}$ and $(1+\epsilon/2)\sigma_{k+1}$. However, aligning with the smaller vectors in this range rather than the larger vectors can incur a cost of at most $O(\epsilon)\sigma_{k+1}^{2}$. Since every column of $\mathbf{Z}$ outside of the first $m$ may incur such a cost, there is a linear accumulation as characterized by Property 2.

Finally, Property 3 captures the intuition that the total error in $\mathbf{Z}$ is bounded by the number of singular values falling in the range $\frac{1}{1+\epsilon/2}\sigma_{k}\leq\sigma_{i}<\sigma_{k}$. This is the total number of singular vectors that aren’t necessarily separated from and can thus be ‘swapped in’ for any of the $(k-m)$ true top vectors with singular value $<(1+\epsilon/2)\sigma_{k+1}$. Property 3 is critical in achieving near optimal Frobenius norm low-rank approximation.

Assume $m\geq 1$. If $m=0$ then Property 1 trivially holds. We will prove the statement for Algorithm 2, since this is the more complex case, and then explain how the proof extends to Algorithm 1.

By Lemma 4 we have with probability $99/100$:

Furthermore, one possible rank $k$ approximation of $p_{1}(\mathbf{A})$ is $p_{1}(\mathbf{A}_{k})$. By the optimality of $p_{1}(\mathbf{A})_{k}$,

The last inequalities follow from setting $q=\Theta(\log(d/\epsilon)/\sqrt{\epsilon})$ and from the fact that $\sigma_{i}\leq\sigma_{k+1}=\alpha$ for all $i\geq k+1$ and thus by property 3 of Lemma 5, $|p_{1}(\sigma_{i})|\leq\frac{\sigma_{k+1}}{2^{q\sqrt{\epsilon/2}-1}}$. Noting that $k\leq d$, we can plug this bound into (4) to get

Applying the Pythagorean theorem and the invariance of the Frobenius norm under rotation gives

Letting $\mathbf{c}_{i}$ be the $i^{\text{th}}$ row of $\mathbf{C}$, expanding out these norms gives

Since $\mathbf{C}$’s columns are orthonormal, its rows all have norms upper bounded by $1$. So $\|\mathbf{c}_{i}\|_{2}^{2}p_{1}(\sigma_{i})^{2}\leq p_{1}(\sigma_{i})^{2}$ for all $i$. So for all $l\leq r$, (6) gives us

Recall that $m$ is the number of singular values with $\sigma_{i}\geq(1+\epsilon/2)\sigma_{k+1}$. By Property 2 of Lemma 5, for all $i\leq m$ we have $\sigma_{i}\leq p_{1}(\sigma_{i})$. This gives, for all $l\leq m$:

Converting these sums back to norms yields $\|\mathbf{\Sigma}_{l}\|_{F}^{2}-\frac{\epsilon\sigma_{k+1}^{2}}{2}\leq\|\mathbf{C}^{T}\mathbf{\Sigma}_{l}\|_{F}^{2}$ and therefore $\|\mathbf{A}_{l}\|_{F}^{2}-\frac{\epsilon\sigma_{k+1}^{2}}{2}\leq\|\mathbf{Y}_{1}\mathbf{Y}_{1}^{T}\mathbf{A}_{l}\|_{F}^{2}$ and

Now $\mathbf{Y}_{1}\mathbf{Y}_{1}^{T}\mathbf{A}_{l}$ is a rank $l$ approximation to $\mathbf{A}$ falling within the column span of $\mathbf{Y}$ and hence within the column span of $\mathbf{Q}$. By Lemma 2, the best rank $l$ Frobenius approximation to $\mathbf{A}$ within $\mathbf{Q}$ is given by $\mathbf{Q}\mathbf{\bar{U}}_{l}(\mathbf{Q}\mathbf{\bar{U}}_{l})^{T}\mathbf{A}$. So we have

For Algorithm 1, we instead choose $p_{1}(x)=(1+\epsilon/2)\sigma_{k+1}\cdot\left(\frac{x}{(1+\epsilon/2)\sigma_{k+1}}\right)^{2q+1}$. For $q=\Theta(\log d/\epsilon)$, this polynomial satisfies the necessary properties: for all $i\geq k+1$, $p_{1}(\sigma_{i})\leq O\left(\frac{\epsilon}{2d^{2}}\sigma_{k+1}^{2}\right)$ and for all $i\leq m$, $\sigma_{i}\leq p_{1}(\sigma_{i})$. Further, up to a rescaling, $p_{1}(\mathbf{A})\mathbf{\Pi}=\mathbf{K}$ so $\mathbf{Y}_{1}$ spans the same space as $\mathbf{K}$. Therefore since Algorithm 1 returns $\mathbf{Z}$ with $\mathbf{Z}_{l}$ equal to the best rank $l$ Frobenius norm approximation to $\mathbf{A}$ within the span of $\mathbf{K}$, for all $l$ we have:

Proof of Property 2

Property 1 and the fact that $\mathcal{E}(\mathbf{Z}_{l},\mathbf{A})$ is always positive immediately gives Property 2 for $l\leq m$. So we need to show that it holds for $m<l\leq k$. Note that if $w$, the number of singular values with $\frac{1}{1+\epsilon/2}\sigma_{k}\leq\sigma_{i}<\sigma_{k}$ is equal to , then $\sigma_{k+1}<\frac{1}{1+\epsilon/2}\sigma_{k}$, so $m=k$ and we are done. So we assume $w\geq 1$ henceforth. Again, we first prove the statement for Algorithm 2 and then explain how the proof extends to the simpler case of Algorithm 1.

Intuitively, Property 1 follows from the guarantee that there is a rank $m$ subspace of $span(\mathbf{K})$ that aligns with $\mathbf{A}$ nearly as well as the space spanned by $\mathbf{A}$’s top $m$ singular vectors. To prove Property 2 we must show that there is also some rank $k$ subspace in $span(\mathbf{K})$ whose components all align nearly as well with $\mathbf{A}$ as $\mathbf{u}_{k}$, the $k^{\text{th}}$ singular vector of $\mathbf{A}$. The existence of such a subspace ensures that $\mathbf{Z}$ performs well, even on singular vectors in the intermediate range $[\sigma_{k},(1+\epsilon/2)\sigma_{k+1}]$.

Let $\mathbf{A}_{inner}$ = $\mathbf{A}_{r\setminus k}-\mathbf{A}_{r\setminus(k+w)}$. $\mathbf{A}_{inner}=\mathbf{U}\mathbf{\Sigma}_{inner}\mathbf{V}^{T}$ where $\mathbf{\Sigma}_{inner}$ contains only the singular values $\sigma_{k+1},\ldots,\sigma_{k+w}$. These are the $w$ intermediate singular values of $\mathbf{A}$ falling in the range $\left[\frac{1}{1+\epsilon/2}\sigma_{k},\sigma_{k}\right)$. Let $\mathbf{A}_{outer}=\mathbf{A}-\mathbf{A}_{inner}=\mathbf{U}\mathbf{\Sigma}_{outer}\mathbf{V}^{T}$. $\mathbf{\Sigma}_{outer}$ contains all large singular values of $\mathbf{A}$ with $\sigma_{i}\geq\sigma_{k}$ and all small singular values with $\sigma_{i}<\frac{1}{1+\epsilon/2}\sigma_{k}$.

We next apply the argument used to prove Property 1 to $p_{2}(\mathbf{A}_{outer})\mathbf{\Pi}$. The $(k+1)^{\text{th}}$ singular value of $\mathbf{A}_{outer}$ is equal to $\sigma_{k+w+1}\leq\frac{1}{1+\epsilon/2}\sigma_{k}=\alpha$. So applying (7) we have for all $l\leq k$,

Plugging (9) and (6.1) into (8) yields that, for any $\mathbf{x}$ in $span(\mathbf{Y}_{2})$, i.e. $span(p_{2}(\mathbf{A})\mathbf{\Pi})$,

So, we have identified a rank $k$ subspace $\mathbf{Y}_{2}$ within our Krylov subspace such that every vector in its span aligns at least as well with $\mathbf{A}$ as $\mathbf{u}_{k}$.

Now, for any $m\leq l\leq k$, consider $\mathcal{E}(\mathbf{Z}_{l},\mathbf{A})$. We know that given $\mathbf{Z}_{l-1}$, we can form a rank $l$ matrix $\mathbf{\overline{Z}}_{l}$ in our Krylov subspace simply by appending a column $\mathbf{x}$ orthogonal to the $l-1$ columns of $\mathbf{Z}_{l-1}$ but falling in the span of $\mathbf{Y}_{2}$. Since $\mathbf{Y}_{2}$ has rank $k$, finding such a column is always possible. Since $\mathbf{Z}_{l}$ is the optimal rank $l$ Frobenius norm approximation to $\mathbf{A}$ falling within our Krylov subspace,

Again, a nearly identical proof applies for Algorithm 1. We just choose $p_{2}(x)=\sigma_{k}\left(\frac{x}{\sigma_{k}}\right)^{2q+1}$. For $q=\Theta(\log d/\epsilon)$ this polynomial satisfies the necessary properties: for all $i\geq k$, $p_{1}(\sigma_{i})\leq O\left(\frac{\epsilon}{2d^{2}}\sigma_{k}^{2}\right)$ and for all $i\leq k$, $\sigma_{i}\leq p_{2}(\sigma_{i})$.

Proof of Property 3

By Properties 1 and 2 we already have, for all $l\leq k$, $\mathcal{E}(\mathbf{Z}_{l},\mathbf{A})\leq\epsilon\sigma_{k+1}^{2}+(l-m)\cdot 3\epsilon\sigma_{k+1}^{2}\leq(1+k-m)\cdot 3\epsilon\cdot\sigma_{k+1}^{2}$. So if $k-m\leq w$ then we immediately have Property 3.

For $l\leq m+w$, then Properties 1 and 2 already give us $\mathcal{E}(\mathbf{Z}_{l},\mathbf{A})\leq\epsilon\sigma_{k+1}^{2}+(l-m)\cdot 3\epsilon\sigma_{k+1}^{2}\leq(w+1)\cdot 3\epsilon\cdot\sigma_{k+1}^{2}$. So consider $m+w\leq l\leq k$. Given $\mathbf{Z}_{m}$, to form a rank $l$ matrix $\mathbf{\overline{Z}}_{l}$ in our Krylov subspace we need to append $l-m$ orthonormal columns. We can choose $\min\{k-w-m,l-m\}$ columns, $\mathbf{X}_{1}$, from the $k-w$ dimensional subspace within $span(\mathbf{Y}_{2})$ that is entirely contained in $span(\mathbf{Y}_{outer})$. If necessary (i.e. $k-w-m\leq l-m$), We can then choose the remaining $l-(k-w)$ columns $\mathbf{X}_{2}$ from the span of $\mathbf{Y}_{2}$.

Similar to our argument when considering a single vector in the span of $\mathbf{Y}_{outer}$, letting $\mathbf{\overline{Y}}_{outer}=\left(\mathbf{I}-\mathbf{X}_{1}\mathbf{X}_{1}^{T}\right)\mathbf{Y}_{outer}$, we have by (10):

By applying (6.1) directly to each column of $\mathbf{X}_{2}$ we also have:

Assume that $\min\{k-w-m,l-m\}=k-w-m$. Similar calculations show the same result when $\min\{k-w-m,l-m\}=l-m$. We can use the above two bounds to obtain:

2 Error Bounds for Simultaneous Iteration and Block Krylov Iteration

With Lemma 9 in place, we can easily prove that Simultaneous Iteration and Block Krylov Iteration both achieve the low-rank approximation and PCA guarantees (1), (2), and (3).

With probability $99/100$, Algorithms 1 and 2 return $\mathbf{Z}$ satisfying (2):

Let $m$ be the number of singular values with $\sigma_{i}\geq(1+\epsilon/2)\sigma_{k+1}$. If $m=0$ then we are done since any $\mathbf{Z}$ will satisfy $\|\mathbf{A}-\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}\|_{2}\leq\|\mathbf{A}\|_{2}=\sigma_{1}\leq(1+\epsilon/2)\sigma_{k+1}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|_{2}$. Otherwise, by Property 1 of Lemma 9,

Additive error in Frobenius norm directly translates to additive spectral norm error. Specifically, applying Theorem 3.4 of , which we also prove as Lemma 15 in Appendix A,

Finally, $\mathbf{Z}_{m}\mathbf{Z}_{m}^{T}\mathbf{A}=\mathbf{Z}\mathbf{Z}_{m}^{T}\mathbf{A}$ and so by Lemma 3 we have $\|\mathbf{A}-\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}\|_{2}^{2}\leq\|\mathbf{A}-\mathbf{Z}_{m}\mathbf{Z}_{m}^{T}\mathbf{A}\|_{2}^{2}$, which combines with (13) to give the result. ∎

With probability $99/100$, Algorithms 1 and 2 return $\mathbf{Z}$ satisfying (1):

$w$ is defined as the number of singular values with $\frac{1}{1+\epsilon/2}\sigma_{k}\leq\sigma_{i}<\sigma_{k}$. So $\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}\geq w\cdot\left(\frac{1}{1+\epsilon/2}\sigma_{k}\right)^{2}$. Plugging into (6.2) we have:

Adjusting constants on the $\epsilon$ gives us the result. ∎

With probability $99/100$, Algorithms 1 and 2 return $\mathbf{Z}$ satisfying (3):

First note that $\mathbf{z}_{i}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{z}_{i}\leq\mathbf{u}_{i}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{u}_{i}$. This is because $\mathbf{z}_{i}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{z}_{i}=\mathbf{z}_{i}^{T}\mathbf{Q}\mathbf{Q}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{Q}\mathbf{Q}^{T}\mathbf{z}_{i}=\sigma_{i}(\mathbf{Q}\mathbf{Q}^{T}\mathbf{A})^{2}$ by our choice of $\mathbf{z}_{i}$. $\sigma_{i}(\mathbf{Q}\mathbf{Q}^{T}\mathbf{A})^{2}\leq\sigma_{i}(\mathbf{A})^{2}$ since applying a projection to $\mathbf{A}$ will decrease each of its singular values (which follows for example from the Courant-Fischer min-max principle). Then by Property 2 of Lemma 9 we have, for all $i\leq k$,

$\sigma_{i}^{2}=\mathbf{u}_{i}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{u}_{i}$, so simply adjusting constants on $\epsilon$ gives the result. ∎

Improved Convergence With Spectral Decay

In order to avoid inverse dependence on the potentially small singular value gap $\frac{\sigma_{k}}{\sigma_{k+1}}-1$, the number of Block Krylov iterations inherently depends on $1/\sqrt{\epsilon}$. This ensures that our matrix polynomial sufficiently separates small singular values from larger ones. However, when $\sigma_{k}>(1+\epsilon)\sigma_{k+1}$ we can actually use $q=\Theta\left(\log(d/\epsilon)/\sqrt{\min\{1,\frac{\sigma_{k}}{\sigma_{k+1}}-1\}}\right)$ iterations, which is sufficient for separating the top $k$ singular values significantly from the lower values. Specifically, if we set $\alpha=\sigma_{k+1}$ and $\gamma=\frac{\sigma_{k}}{\sigma_{k+1}}-1$, we know that with $q=\Theta\left(\log(d/\epsilon)/\sqrt{\min\{1,\frac{\sigma_{k}}{\sigma_{k+1}}-1\}}\right)$, (5) still holds. We can then just follow the proof of Lemma 9 and show that Property 1 holds for all $l\leq k$ (not just for $l\leq m$ as originally proven). This gives Property 2 and Property 3 trivially.

Further, for $p\geq k$, the exact same analysis shows that $q=\Theta\left(\log(d/\epsilon)/\sqrt{\min\{1,\frac{\sigma_{k}}{\sigma_{p+1}}-1\}}\right)$ suffices. When $\mathbf{A}$’s spectrum decays rapidly, so $\sigma_{p+1}\leq c\cdot\sigma_{k}$ for some constant $c<1$ and some $p$ not much larger than $k$, we can obtain significantly faster runtimes. Our $\epsilon$ dependence becomes logarithmic, rather than polynomial:

With probability $99/100$, for any $p\geq k$, Algorithm 1 or 2 initialized with $\mathbf{\Pi}\sim\mathcal{N}(0,1)^{d\times p}$ returns $\mathbf{Z}$ satisfying guarantees (1), (2), and (3) as long as we set $q=\Theta\left(\log(d/\epsilon)/\left(\min\{1,\frac{\sigma_{k}}{\sigma_{p+1}}-1\}\right)\right)$ or $\Theta\left(\log(d/\epsilon)/\sqrt{\min\{1,\frac{\sigma_{k}}{\sigma_{p+1}}-1\}}\right)$, respectively.

This theorem may prove especially useful in practice because, on many architectures, multiplying a large $\mathbf{A}$ by $2k$ or even $10k$ vectors is not much more expensive than multiplying by $k$ vectors. Additionally, it should still be possible to perform all steps for post-processing $\mathbf{K}$ in memory, again limiting additional runtime costs due to its larger size.

Finally, we note that while Theorem 13 is more reminiscent of classical gap-dependent bounds, it still takes substantial advantage of the fact that we’re looking for nearly optimal low-rank approximations and principal components instead of attempting to converge precisely to $\mathbf{A}$’s true singular vectors. This allows the result to avoid dependence on the gap between adjacent singular values, instead varying only with $\frac{\sigma_{k}}{\sigma_{p+1}}$, which should be much larger.

Experiments

We close with several experimental results. A variety of empirical papers, not to mention widespread adoption, already justify the use of randomized SVD algorithms. Prior work focuses in particular on benchmarking Simultaneous Iteration and, due to its improved accuracy over sketch-and-solve approaches, this algorithm is popular in practice . As such, we focus on demonstrating that for many data problems Block Krylov Iteration can offer significantly better convergence.

We implement both algorithms in MATLAB using Gaussian random starting matrices with exactly $k$ columns. We explicitly compute $\mathbf{K}$ for both algorithms, as described in Section 5, and use reorthonormalization at each iteration to improve stability . We test the algorithms with varying iteration count $q$ on three common datasets, SNAP/amazon0302 , SNAP/email-Enron , and 20 Newsgroups , computing column principal components in all cases. We plot error vs. iteration count for metrics (1), (2), and (3) in Figure 3. For per vector error (3), we plot the maximum deviation amongst all top $k$ approximate principal components (relative to $\sigma_{k+1}$).

Unsurprisingly, both algorithms obtain very accurate Frobenius norm error, $\|\mathbf{A}-\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}\|_{F}/\|\mathbf{A}-\mathbf{A}_{k}\|_{F}$, with very few iterations. This is our intuitively weakest guarantee and, in the presence of a heavy singular value tail, both iterative algorithms will outperform the worst case analysis.

On the other hand, for spectral norm low-rank approximation and per vector error, we confirm that Block Krylov Iteration converges much more rapidly than Simultaneous Iteration, as predicted by our theoretical analysis. It it often possible to achieve nearly optimal error with $<8$ iterations where as getting to within say $1\%$ error with Simultaneous Iteration can take much longer.

The final plot in Figure 3 shows error verses runtime for the $11269\times 15088$ dimensional 20 Newsgroups dataset. We averaged over 7 trials and ran the experiments on a commodity laptop with 16GB of memory. As predicted, because its additional memory overhead and post-processing costs are small compared to the cost of the large matrix multiplication required for each iteration, Block Krylov Iteration outperforms Simultaneous Iteration for small $\epsilon$.

More generally, these results justify the importance of convergence bounds that are independent of singular value gaps. Our analysis in Section 7 predicts that, once $\epsilon$ is small in comparison to the gap $\frac{\sigma_{k}}{\sigma_{k+1}}-1$, we should see much more rapid convergence since $q$ will depend on $\log(1/\epsilon)$ instead of $1/\epsilon$. However, for Simultaneous Iteration, we do not see this behavior with SNAP/amazon0302 and it only just begins to emerge for 20 Newsgroups.

While all three datasets have rapid singular value decay, a careful look confirms that their singular value gaps are actually quite small! For example, $\frac{\sigma_{k}}{\sigma_{k+1}}-1$ is .004 for SNAP/amazon0302 and .011 for 20 Newsgroups, in comparison to .042 for SNAP/email-Enron. Accordingly, the frequent claim that singular value gaps can be taken as constant is insufficient, even for small $\epsilon$.

We thank David Woodruff, Aaron Sidford, Richard Peng and Jon Kelner for several valuable conversations. Additionally, Michael Cohen was very helpful in discussing many details of this project, including the ultimate form of Lemma 9. This work was partially supported by NSF Graduate Research Fellowship Grant No. 1122374, AFOSR grant FA9550-13-1-0042, DARPA grant FA8650-11-C-7192, and the NSF Center for Science of Information.

References

Appendix A Appendix

We first give a deterministic Lemma, from which the main approximation result follows.

We follow . Apply Lemma 14 with $\mathbf{S}=\mathbf{\Pi}$. With probability $1$, $\mathbf{V}_{k}^{T}\mathbf{S}$ has full rank. So, to show the result we need to show that $\|\left(\mathbf{A}-\mathbf{A}_{k}\right)\mathbf{S}\left(\mathbf{V}_{k}^{T}\mathbf{S}\right)^{+}\|_{F}^{2}\leq c\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}$ for some fixed $c$. For any two matrices $\mathbf{M}$ and $\mathbf{N}$, $\|\mathbf{MN}\|_{F}\leq\|\mathbf{M}\|_{F}\|\mathbf{N}\|_{2}$. This property is known as spectral submultiplicativity. Noting that $\|\mathbf{U}_{r\setminus k}\mathbf{\Sigma}_{r\setminus k}\|_{F}^{2}=\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}$ and applying submultiplicativity,

By the rotational invariance of the Gaussian distribution, since the rows of $\mathbf{V}^{T}$ are orthonormal, the entries of $\mathbf{V}_{k}^{T}\mathbf{S}$ and $\mathbf{V}^{T}_{r\setminus k}\mathbf{S}$ are independent Gaussians. By standard Gaussian matrix concentration results (Fact 6 of , also in ), with probability at least $99/100$, $\|\mathbf{V}^{T}_{r\setminus k}\mathbf{S}\|_{2}^{2}\leq c_{1}\cdot\max\{k,r-k\}\leq c_{1}\dot{d}$ and $\|\left(\mathbf{V}_{k}^{T}\mathbf{S}\right)^{+}\|_{2}^{2}\leq c_{2}k$ for some fixed constants $c_{1},c_{2}$. So,

for some fixed $c$, yielding the result. Note that we choose probability $99/100$ for simplicity – we can obtain a result with higher probability by simply allowing for a higher constant $c$, which in our applications of Lemma 4 will only factor into logarithmic terms. ∎

Chebyshev Polynomials

Given a specified value $\alpha>0$, gap $\gamma\in(0,1]$, and $q\geq 1$, there exists a degree $q$ polynomial $p(x)$ such that:

$p(\left(1+\gamma)\alpha\right)=(1+\gamma)\alpha$

$p(x)\geq x$ for all $x\geq(1+\gamma)\alpha$

$|p(x)|\leq\frac{\alpha}{2^{q\sqrt{\gamma}-1}}$ for all $x\in[0,\alpha]$

Furthermore, when $q$ is odd, the polynomial only contains odd powered monomials.

The required polynomial can be constructed using a standard Chebyshev polynomial of degree $q$, $T_{q}(x)$, which is defined by the three term recurrence:

Each Chebyshev polynomial satisfies the well known property that $T_{q}(x)\leq 1$ for all $x\in$ and, for $x>1$, we can write the polynomials in closed form :

which is clearly of degree $q$ and well defined since, referring to (15), $T_{q}(x)>0$ for all $x>1$. Now,

so $p(x)$ satisfies property 1. With property 1 in place, to prove that $p(x)$ satisfies property 2, it suffices to show that $p^{\prime}(x)\geq 1$ for all $x\geq(1+\gamma)\alpha$. By chain rule,

Thus, it suffices to prove that, for all $x\geq(1+\gamma)$,

We do this by showing that $(1+\gamma)T^{\prime}_{q}(1+\gamma)\geq T_{q}(1+\gamma)$ and then claim that $T^{\prime\prime}_{q}(x)\geq 0$ for all $x>(1+\gamma)$, so (17) holds for $x>(1+\gamma)$ as well. A standard form for the derivative of the Chebyshev polynomial is