Subspace Iteration Randomization and Singular Value Problems

Introduction

Randomized algorithms have established themselves as some of the most competitive methods for rapid low-rank matrix approximation, which is vital in many areas of scientific computing, including principal component analysis and face recognition , large scale data compression and fast approximate algorithms for PDEs and integral equations . In this paper, we consider randomized algorithms for low-rank approximations and singular value approximations within the subspace iteration framework, leading to results that simultaneously retain the reliability of randomized algorithms and the typical faster convergence of subspace iteration methods.

Given any $m\times n$ matrix $A$ with $m\geq n$, its singular value decomposition (SVD) is described by the equation

where $U$ is an $m\times n$ column orthogonal matrix; $V$ is an $n\times n$ orthogonal matrix; and $\Sigma=\mathop{\operator@font diag}\nolimits(\sigma_{1},\cdots,\sigma_{n})$ with $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{n}\geq 0$. Writing $U$ and $V$ in terms of their columns,

then $u_{j}$ and $v_{j}$ are the left and right singular vectors corresponding to $\sigma_{j}$, the $j$-th largest singular value of $A$. For any $1\leq k\leq n$, we let

be the (rank-$k$) truncated SVD of $A$. The matrix $A_{k}$ is unique only if $\sigma_{k+1}<\sigma_{k}$. The assumption that $m\geq n>\max{k,2}$ will be maintained throughout this paper for ease of exposition. Our results still hold for $m<n$ by applying all the algorithms on $A^{T}$. Similarly, all our main results are derived under the assumption that ${\bf rank}(A)=n$. But they remain unchanged even if ${\bf rank}(A)<n$, and hence remain valid by a continuity argument. All our analysis is done without consideration of round-off errors, and thus need not hold exactly true in finite precision, especially when the user tolerances for the low-rank approximation are close to machine precision levels. Additionally, we assume throughout this paper that all matrices are real. In general, $A_{k}$ is an ideal rank-$k$ approximation to $A$, due to the following celebrated property of the SVD:

While there are results similar to Theorem 1 for all unitarily invariant matrix norms, our work on low-rank matrix approximation bounds will only focus on the two most popular of such norms: the 2-norm and the Frobenius norm.

Theorem 1 states that the truncated SVD provides a rank-$k$ approximation to $A$ with the smallest possible 2-norm error and Frobenius-norm error. In the 2-norm, any rank-$k$ approximation will result in an error no less than $\sigma_{k+1}$, and in the Frobenius-norm, any rank-$k$ approximation will result in an error no less than $\sqrt{\sum_{j=k+1}^{n}\sigma_{j}^{2}}$. Additionally, the singular values of $A_{k}$ are exactly the first $k$ singular values of $A$, and the singular vectors of $A_{k}$ are the corresponding singular vectors of $A$. Note, however, that while the solution to problem (3) must be $A_{k}$, solutions to problem (2) are not unique and include, for example, the rank-$k$ matrix $B$ defined below for any $0\leq\theta\leq 1$:

This subtle distinction between the 2-norm and Frobenius norm will later on become very important in our analysis of randomized algorithms (see Remark 3.4.) In Theorem 8 we prove an interesting result related to Theorem 1 for rank-$k$ approximations that only solve problems (2) and (3) approximately.

To compute a truncated SVD of a general $m\times n$ matrix $A$, one of the most straightforward techniques is to compute the full SVD and truncate it, with a standard linear algebra software package like the LAPACK . This procedure is stable and accurate, but it requires $O(mn^{2})$ floating point operations, or flops. This is prohibitively expensive for applications such as data mining, where the matrices involved are typically sparse with huge dimensions. In other practical applications involving the truncated SVD, often the very objective of computing a rank-$k$ approximation is to avoid excessive computation on $A$. Hence it is desirable to have schemes that can compute a rank-$k$ approximation more efficiently. Depending on the reliability requirements, a good rank-$k$ approximation can be a matrix that is accurate to within a constant factor from the optimal, such as a rank-revealing factorization (more below), or it can be a matrix that closely approximates the truncated SVD itself.

Many approaches have been taken in the literature for computing low-rank approximations, including rank-revealing decompositions based on the QR, LU, or two-sided orthogonal (aka UTV) factorizations . Recently, there has been an explosion of randomized algorithms for computing low-rank approximations . There is also software package available for computing interpolative decompositions, a form of low-rank approximation, and for computing the PCA, with randomized sampling . These algorithms are attractive for two main reasons: they have been shown to be surprisingly efficient computationally; and like subspace methods, the main operations involved in many randomized algorithms can be optimized for peak machine performance on modern architectures. For a detailed analysis of randomized algorithms and an extended reference list, see ; for a survey of randomized algorithms in data analysis, see .

The subspace iteration is a classical approach for computing singular values. There is extensive convergence analysis on subspace iteration methods and a large literature on accelerated subspace iteration methods . In general, it is well-suited for fast computations on modern computers because its main computations are in terms of matrix-matrix products and QR factorizations that have been highly optimized for maximum efficiency on modern serial and parallel architectures . There are two well-known weaknesses of subspace iteration, however, that limit its practical use. On one hand, subspace iteration typically requires very good separation between the wanted and unwanted singular values for good convergence. On the other hand, good convergence also often critically depends on the choice of a good start matrix .

Another classical class of approximation methods for computing an approximate SVD are the Krylov subspace methods, such as the Lanczos algorithm (see, for example .) The computational cost of these methods depends heavily on several factors, including the start vector, properties of the input matrix and the need to stabilize the algorithm. One of the most important part of the Krylov subspace methods, however, is the need to do a matrix-vector product at each iteration. In contrast to matrix-matrix products, matrix-vector products perform very poorly on modern architectures due to the limited data reuse involved in such operations, In fact, one focus of Krylov subspace research is on effective avoidance of matrix-vector operations in Krylov subspace methods (see, for example .)

This work focuses on the surprisingly strong performance of randomized algorithms in delivering highly accurate low-rank approximations and singular values. To illustrate, we introduce Algorithm 1.1, one of the basic randomized algorithms (see .)

Compute an orthogonal column basis $Q$ for $Y$.

Compute $B_{k}$, the rank-$k$ truncated SVD of $B$.

Throughout this paper, a random matrix, such as $\Omega$ in Algorithm 1.1, is a standard Gaussian matrix, i.e., its entries are independent standard normal variables of zero mean and standard deviation $1$.

While other random matrices might work equally well, the choice of the Gaussian matrix provides two unique advantages: First, the distribution of a standard Gaussian matrix is rotationally invariant: If $V$ is an orthonormal matrix, then $V^{T}\Omega$ is itself a standard Gaussian matrix with the same statistical properties as $\Omega$ . Second, our analysis is much simplified by the vast literature on the singular value probability density functions of the Gaussian matrix.

While Algorithm 1.1 looks deceptively simple, its analysis is long, arduous, and involves very strong doses of statistics . The following theorem establishes an error bound on the accuracy of $QQ^{T}A$ as a low-rank approximation to $A$. There are similar results in the Frobenius norm.

(Halko, Martinsson, Tropp [35, Corollary 10.9]) The column-orthonormal matrix $Q$ produced by Step 3 in Algorithm 1.1 satisfies

with failure probability at most $6e^{-p}$.

Comparing Theorem 2 with Theorem 1, it is clear that Algorithm 1.1 could provide a very good low rank approximation to $A$ with probability at least $1-6e^{-p}$, despite its simple operations, provided that $\sigma_{k+1}\ll\|A\|_{2}$. While algorithms differ in their algorithm design, efficiency, and domain applicability, they typically share the same advantages of computational efficiency and approximation accuracy.

Algorithm 1.1 is the combination of Stages A and B of the Proto Algorithm in , where the truncated SVD is considered separately from low-rank approximation. In Section 2.3 we will discuss the pros and cons of SVD truncation vs. no truncation. Algorithm 1.1 is a special case of the randomized subspace iteration method (see Algorithm 2.2), for which Halko, Martinsson, Tropp have developed similar results.

However, while the upper bound in Theorem 2 can be very satisfactory for many applications, there may be situations where singular value approximations are also desirable. In addition, it is well-known that in practical computations randomized algorithms often far outperform their error bounds , whereas the results in do not suggest convergence of the computed rank-$k$ approximation to the truncated SVD in either Algorithm 1.1 or the more general randomized subspace iteration method.

Our entire work is based on novel analysis of the subspace iteration method, and we consider randomized algorithms within the subspace iteration framework. This allows us to take advantage of existing theories and technical machinery in both fields.

Current analysis on randomized algorithms focuses on the errors in the approximation of $A$ by a low rank matrix, whereas classical analysis on subspace iteration methods focuses on the accuracy in the approximate singular values. Our analysis allows us to obtain both kinds of results for both of these methods, leading to the stronger rank-revealing approximations. In terms of randomized algorithms, our matrix approximation bounds are in general tighter and can be drastically better than existing ones; in terms of singular values, our relative convergence lower bounds can be interpreted as simultaneously convergence error bounds and rank-revealing lower bounds.

Our analysis has lead us to some interesting conclusions, all with high probability (more precise statements are in Sections 5 through 7):

The leading $k$ singular values computed by randomized algorithms are at least a good fraction of the true ones, regardless of how the singular values are distributed, and they converge quickly to the true singular values in case of rapid singular value decay. In particular, this result implies that randomized algorithms can also be used as efficient and reliable condition number estimators.

The above results, together with the fact that randomized algorithms compute low-rank approximations up to a dimension dependent constant factor from optimal, mean that these low-rank approximations are in fact rank-revealing factorizations. In addition, for rapidly decaying singular values, these approximations can be as accurate as a truncated SVD.

The subspace iteration method in general and the power method in particular is still slowly convergent without over-sampling in the start matrix. We present an alternative choice of the start matrix based on our analysis, and demonstrate its competitiveness.

The rest of this paper is organized as follows: In Section 2 we discuss subspace iteration methods and their randomized versions in more detail; in Section 3 we list a number of preliminary as well as key results needed for later analysis; in Section 4 we derive deterministic lower bounds on singular values and upper bounds on low-rank approximations; in Section 5 we provide both average case and large deviation bounds on singular values and low-rank approximations; in Section 6 we compare these approximations with other rank-revealing factorizations; in Section 7 we discuss how randomized algorithms can be used as efficient and reliable condition number estimators; in Section 8 we present supporting numerical experimental results; and in Section 9 we draw some conclusions and point out possible directions for future research.

Much of our analysis has its origin in the analysis of subspace iteration and randomized algorithms . It relies both on linear algebra tools as well as statistical analysis to do some of the needed heavy lifting to reach our conclusions. To limit the length of this paper, we have put the more detailed parts of the analysis as well as some additional numerical experimental results in the Supplemental Material, which is accessible at SIAM’s on-line portal.

Algorithms

In this section, we present the main algorithms that are discussed in the rest of this paper. We also discuss subtle differences between our presentation of randomized algorithms and that in .

We start with the classical subspace iteration method for computing the largest few singular values of a given matrix.

Compute $Y=\left(AA^{T}\right)^{q}A\,\Omega$.

Compute an orthogonal column basis $Q$ for $Y$.

Compute $B_{k}$, the rank-$k$ truncated SVD of $B$.

Given the availability of Lanczos-type algorithms for the singular value computations, the classical subspace iteration method is not widely used in practice except when $k\ll n$. We present it here for later comparisons with its randomized version. We ignore the vast literature of accelerated subspace iteration methods (see, for example ) in this paper since our main goal here is to analyze the convergence behavior of subspace iteration method with and without randomized start matrix $\Omega$.

We have presented Algorithm 2.1 in an over-simplified form above to convey the basic ideas involved. In practice, the computation of $Y$ would be prone to round-off errors. For better numerical accuracy, Algorithm A.1 in the Appendix should be used numerically to compute the $Q$ matrix in Algorithm 2.1. In practical computations, however, Algorithm A.1 is often performed once every few iterations, to balance efficiency and numerical stability (see Saad .) In the rest of Section 2, any QR factorization of the matrix $Y=\left(AA^{T}\right)^{q}A\Omega$ should be computed numerically through periodic use of Algorithm A.1.

While there is little direct analysis of subspace iteration methods for singular values (Algorithm 2.1) in the literature, one can generalize results of subspace iteration methods for symmetric matrices to the singular value case in a straightforward fashion. The symmetric matrix version of Theorem 3 can be found in .

(Bathe and Wilson) Assume that Algorithm 2.1 converges as $q\rightarrow\infty$. Then

2 Randomized Algorithms

In order to enhance the convergence of Algorithm 2.1 in the absence of any useful information about the leading singular vectors, a sensible approach is to replace the deterministic start matrix with a random one, leading to

Compute a rank-$k$ approximation with Algorithm 2.1.

Since Algorithm 2.2 is the special case of Algorithm 2.1 with $\Omega$ being chosen as random, all our results for Algorithm 2.1 equally hold for Algorithm 2.2.

The only difference between Algorithm 2.1 and Algorithm 2.2 is in the choice of $\Omega$, yet this difference will lead to drastically different convergence behavior. One of the main purposes of this paper is to show that the slow or non-convergence of Algorithm 2.1 due to bad choice of $\Omega$ vanishes with near certainty in Algorithm 2.2. In particular, a single iteration ($q=0$ in Algorithm 2.2) in the randomized subspace iteration method is often sufficient to return good enough singular values and low-rank approximations (Section 5).

Our analysis of deterministic and randomized subspace iteration method was in large part motivated by the analysis and discussion of randomized algorithms in . We have chosen to present the algorithms in Section 2 in forms that are not identical to those in for ease of stating our results in Sections 4 through 8. Versions of Algorithm 2.2 have also appeared in for solving large-scale discrete inverse problems.

3 To Truncate or not to Truncate

The randomized algorithms in Section 2 are presented in a slight different form than those in . One key difference is in the step of SVD truncation, which is considered an optional postprocessing step there. In this section, we discuss the pros and cons of SVD truncation. We start with the following simple lemma, versions of which appear in .

Lemma 4 makes it obvious that any SVD truncation of $Q^{T}A$ will only result in a less accurate approximation in the 2-norm and Frobenius norm. This is strong motivation for no SVD truncation. The SVD truncation of $Q^{T}A$ also involves the computation of the SVD of $Q^{T}A$ in some form, which also results in extra computation.

Setup

In this section we build some of the technical machinery needed for our heavy analysis later on. We start by reciting two well-known results in matrix analysis, and then develop a number of theoretical tools that outline our approach in the low-rank approximation analysis. Some of these results may be of interest in their own right. For any matrix $X$, we use $\sigma_{j}(X)$ to denote its $j$-th largest singular value.

The Cauchy interlacing theorem shows the limitations of any approximation with an orthogonal projection.

(Golub and van Loan [30, p. 411]) Let $A$ be an $m\times n$ matrix and $Q$ be a matrix with orthonormal columns. Then $\sigma_{j}(A)\geq\sigma_{j}(Q^{T}A)$ for $1\leq j\leq\min(m,n)$.

A direct consequence of Theorem 5 is that $\sigma_{j}(A)\geq\sigma_{j}(\widehat{A})$, where $\widehat{A}$ is any submatrix of $A$.

Weyl’s monotonicity theorem relates singular values of matrices $X$ and $Y$ to those of $X+Y$.

(Weyl’s monotonicity theorem [43, Thm. 3.3.16]) Let $X$ and $Y$ be $m\times n$ matrices with $m\geq n$. Then

The Hoffman-Wielandt theorem bounds the errors in the differences between the singular values of $X$ and those of $Y$ in terms of $\|X-Y\|_{F}$.

(Hoffman and Wielandt ) Let $X$ and $Y$ be $m\times n$ matrices with $m\geq n$. Then

Below we develop a number of theoretical results that will form the basis for our later analysis on low-rank approximations. Theorem 8 below is of potentially broad independent interest. Let $B$ be a rank-$k$ approximation to $A$. Theorem 8 below relates the approximation error in the Frobenius norm to that in the 2-norm as well as the approximation errors in the leading $k$ singular values. It will be called the Reverse Eckart and Young Theorem due to its complimentary nature with Theorem 1 in the Frobenius norm.

(Reverse Eckart and Young) Given any $m\times n$ matrix $A$, and let $B$ be a matrix with rank at most $k$ such that

when $\eta$ is larger than or close to $\sigma_{k+1}$. On the other hand, if $\eta\ll\sigma_{k+1}$, then equation (6) simplifies to

where the last ratio can be much smaller than $\eta$, implying a much better rank-$k$ approximation in $B$. Similar comments apply to equation (5). This interesting feature of Theorem 8 is one of the reasons why our eventual 2-norm and Frobenius norm upper bounds are much better than those in Theorem 2 in the event that $\eta\ll\sigma_{k+1}$. This also has made our proofs in Appendix B somewhat involved in places.

Equation (7) asserts that a small $\eta$ in equation (5) necessarily means good approximations to all the $k$ leading singular values of $A$. In particular, $\eta=0$ means the leading $k$ singular values of $A$ and $B$ must be the same. However, our singular value analysis will not be based on Equation (7), as our approach in Section 4 provides us with much better results.

Proof of Theorem 8: Write $A=\left(A-B\right)+B$. It follows from Theorem 6 that for any $1\leq i\leq n-k$:

since $B$ is a rank-$k$ matrix. It follows that

Plugging this into equation (5) yields (6).

As to equation (7), we observe that the $(k+1)-st$ through the last singular values of $B$ are all zero, given that $B$ has rank $k$. Hence the result trivially follows from Theorem 7,

Our next theorem is a generalization of Theorem 1.

Problem (9) in Theorem 9 is a type of restricted SVD problem. Oddly enough, this problem becomes much harder to solve for the 2-norm. In fact, $B_{k}$ might not even be the solution to the corresponding restricted SVD problem in 2-norm. Combining Theorems 8 and 9, we obtain

Our low-rank approximation analysis in the $2$-norm will be based on equation (11). While this is sufficient, it also makes our $2$-norm results perhaps weaker than they should be due to the mixture of the $2$-norm and the Frobenius norm.

By Theorem 1, $A_{k}$ is the best Frobenius norm approximation to $A$, whereas by Theorem 9 $QB_{k}$ is the best restricted Frobenius norm approximation to $A$. This leads to the following interesting consequence

Thus we can expect $QB_{k}$ to also be an excellent rank-$k$ approximation to $A$ as long as $Q$ points to the principle singular vector directions.

Result (9) is now an immediate consequence of Theorem 1. To prove (10), we observe that

The third term in the last equation is zero because $A_{k}\left(A-A_{k}\right)^{T}=0.$ Combining this last relation with equation (12) gives us relation (10). Q.E.D.

Deterministic Analysis

In this section we perform deterministic convergence analysis on Algorithm 2.1. Theorem 12 is a relative convergence lower bound, and Theorem 13 is an upper bound on the matrix approximation error. Both appear to be new for subspace iteration. Our approach, while quite novel, was motivated in part by the analysis of subspace iteration methods by Saad and randomized algorithms in . Since Algorithm 1.1 is a special case of Algorithm 2.2 with $q=0$, which in turn is a special case of Algorithm 2.1 with an initial random matrix, our analysis applies to them as well and will form the basis for additional probabilistic analysis in Section 5.

We begin by noticing that the output $QB_{k}$ in Algorithm 2.1 is also the rank-$k$ truncated SVD of the matrix $QQ^{T}A$, due to the fact that $Q$ is column orthonormal. In fact, columns of $Q$ are nothing but an orthonormal basis for the column space of matrix $\left(AA^{T}\right)^{q}A\,\Omega$. This is the reason why Algorithm 2.1 is called subspace iteration. Lemma 10 below shows how to obtain alternative orthonormal bases for the same column space. We omit the proof.

The matrix ${\displaystyle\widehat{\Omega}_{1}}$ has at least as many columns as rows. Assume it is of full row rank so that its pseudo-inverse satisfies

Below we present a special choice of $X$ that will reveal the manner in which convergence to singular values and low-rank approximations takes place. Ideally, such an $X$ would orient the first $k$ columns of {\displaystyle U\left(\begin{array}[]{c}\left(\begin{array}[]{cc}\Sigma_{1}&\cr&\Sigma_{2}\end{array}\right)^{2q+1}\widehat{\Omega}_{1}\cr\cr\Sigma_{3}^{2q+1}\widehat{\Omega}_{2}\end{array}\right)X} in the directions of the leading $k$ singular vectors in $U$. We choose

By equation (16), the QR factorization of $\left(AA^{T}\right)^{q}A\Omega X$ can now be written in the following $3\times 3$ partition:

We will use this representation to derive convergence upper bounds for singular value and rank-k approximations. In particular, we will make use of the fact that the above QR factorization also embeds another one

We are now ready to derive a lower bound on $\sigma_{k}(B_{k})$.

Let $H_{1}$ be defined in equation (16), and assume that the matrix $\widehat{\Omega}_{1}$ has full row rank, then the matrix $B_{k}$ computed in Algorithm 2.1 must satisfy

It might seem more intuitive in equation (15) to choose $X=\left(X_{1}\quad X_{2}\right)$ where $X_{1}$ solves the following least squares problem

Our choice of $X$ seems as effective and allows simpler analysis.

Proof of Lemma 11: We note by Lemma 10 that

From equations (20) and (18), we see that the matrix

is simply a submatrix of the middle matrix on the right hand side of equation (20). By Remark 3.1, it follows immediately that

Combining these two relations, and together with the fact that ${\displaystyle\left\|\widehat{R}_{11}^{T}\right\|_{2}=\sqrt{1+\left\|H_{1}\right\|_{2}^{2}}}$, we obtain (19). Q.E.D.

2 Deterministic Bounds

In this section we develop the analysis in Section 4.1 into deterministic lower bounds for singular values and upper bounds for rank-k approximations.

Since the interlacing theorem 5 asserts an upper bound $\sigma_{k}\left(B_{k}\right)\leq\sigma_{k}$, equation (19) provides a nice lower bound on $\sigma_{k}\left(B_{k}\right)$. These bounds mean that $\sigma_{k}\left(B_{k}\right)$ is a good approximation to $\sigma_{k}$ as long as $\left\|H_{1}\right\|_{2}$ is small. This consideration is formalized in the theorem below.

Proof of Theorem 12: By the definition of the matrix $H_{1}$ in equation (16), it is straightforward to get

This, together with lower bound (19), gives the result in Theorem 12 for $j=k$. To prove Theorem 12 for any $1\leq j<k$, we observe that since $\sigma_{j}\left(B_{k}\right)=\sigma_{j}\left(B_{j}\right)$, all that is needed is to repeat all previous arguments for a rank $j$ truncated SVD. Q.E.D.

Now we consider rank-k approximation upper bounds. Toward this end and considering Theorem 9, we would like to start with an upper bound on $\|\left(I-QQ^{T}\right)A_{k}\|_{F}$. By Lemma 10 and equation (17), we have

Since $A_{k}=U\mathop{\operator@font diag}\nolimits\left(\Sigma_{1},0,0\right)V^{T}$, and since \widehat{Q}_{1}={\displaystyle U\left(\begin{array}[]{c}I\cr 0\cr H_{1}\end{array}\right)\widehat{R}_{11}^{-1}} according to equation (18), the above right hand side becomes

where we have used the fact that (see (18))

However, when $\Omega$ is taken to be Gaussian, only the bounds in Theorem 13 allow average case analysis for all values of $p$ (See Section 5.2.)

Not surprisingly, the two key factors governing the singular value convergence of Algorithm 2.1 also govern the convergence of the low-rank approximation. Hence Remark 4.2 applies equally well to Theorem 13.

Proof of Theorem 13: We first assume that the matrix $H_{1}$ in equation (16) has full column rank. Rewrite

The last expression is a symmetric positive semi-definite matrix. This allows us to write

By a continuity argument, the last relation remains valid even if $H_{1}$ does not have full column rank.

Due to the special form of $H_{1}$ in equation (16), we can write $H_{1}\Sigma_{1}$ as

Plugging this into equation (34) and dividing both the nomerator and denominator by $\sigma_{1}$,

Comparing this with Theorems 8 and 9 proves Theorem 13. Q.E.D.

Statistical Analysis

This section carries out the needed statistical analysis to reach our approximation error bounds. In Section 5.1 we make a list of the statistical tools used in this analysis; in Section 5.2 we perform average value analysis on our error bounds; and in Section 5.3 we provide large deviation bounds.

The simplest of needed statistical results necessary for our analysis is the following proposition from .

For fix matrices $S,T$ and standard Gaussian matrix $G$, we have

The following large deviation bound for the pseudo-inverse of a Gaussian matrix is also from .

The following theorem provides classical tail bounds for functions of Gaussian matrices. It was taken from [Thm. 4.5.7].

Suppose that $h$ is a real valued Lipschitz function on matrices:

Draw a standard Gaussian matrix $G$. Then

The two propositions below will be used in our average case error bounds analysis, both for singular values and rank-k approximations. Their proofs are lengthy and can be found in the Supplemental Material.

Let $\alpha>0$, $\beta>0$, $\gamma>0$ and $\delta>0$, and let $G$ be an $m\times n$ Gaussian matrix. Then

There are lower and upper bounds similar to Proposition 17 for the pseudo-inverse of a Gaussian, with a significant complication. When $G$ is a square Gaussian matrix, it is non-singular with probability $1$. However, the probability density function for its pseudo-inverse could have a very long tail according to Lemma 15. A similar argument could also be made when $G$ is almost a square matrix. This complication will have important implications for parameter choices in Algorithm 2.2 (see Sections 5.2 and 5.3.) Function $\log(\cdot)$ below is base-$e$.

2 Average Case Error Bounds

This section is devoted to the average case analysis of Algorithm 2.2. This work requires us to study the average case behavior on the upper and lower bounds in Theorems 12 and 13. As observed in Section 2.2, the distribution of a standard Gaussian matrix is rotationally invariant, and hence the matrices $\widehat{\Omega}_{1}$ and $\widehat{\Omega}_{2}$ are themselves independent standard Gaussian matrices. With the tools established in Section 5.1, our analysis here consists mostly of stitching together the right pieces from there.

Let $A=U\Sigma V^{T}$ be the SVD of $A$, and let $QB_{k}$ be a rank-$k$ approximation computed by Algorithm 2.2. Then for $j=1,\cdots,k,$

Theorem 19 strongly suggests that in general some over-sampling in the number of columns can significantly improve convergence in the singular value approximation. This is consistent with the literature and is very significant for practical implementations.

Since $\tau_{j}\leq 1$ for all $j$, Theorem 19 implies that for $p\geq 2$ and for all $j\leq k$,

In other words, Algorithm 2.2 approximates the leading $k$ singular values by a good fraction on average, regardless of how the singular values are distributed, even for $q=0$. This result is surprising and yet valuable. It will have applications in condition number estimation (see Sections 5.3 and 7 for more discussion.)

For matrices with rapidly decaying singular values, convergence could be so rapid that one could even set $q=0$ in some cases (Section 5.3.) This is the basis of the excitement about Algorithm 2.2 in that very little work is typically sufficient to realize an excellent low-rank approximation. The faster the singular values decay, the faster Algorithm 2.2 converges.

Proof of Theorem 19: As in Theorem 12, we will only prove Theorem 19 for $j=k$. All other values of $j$ can be proved by simply citing Theorem 19 for a rank-$j$ SVD truncation. Since $\widehat{\Omega}_{2}$ and $\widehat{\Omega}_{1}$ are independent of each other, we will take expectations over $\widehat{\Omega}_{2}$ and $\widehat{\Omega}_{1}$ in turn, based on Propositions 17 and 18.

Let ${\displaystyle\alpha=\left\|\widehat{\Omega}_{1}^{\dagger}\right\|_{2}\tau_{k}^{2q+1}}$. By Theorem 12 and Proposition 17,

For $p\geq 2$, we further take expectation over $\widehat{\Omega}_{1}$ according to Proposition 18. By equation (50),

To complete the proof, we note that the results for $p=1$ and $p=0$ can be obtained similarly by taking expectation of $\widehat{\Omega}_{1}$ over equation (50) and simplifying. Q.E.D.

It is now time for average case analysis of low-rank matrix approximations. Again, we base our arguments on Propositions 17 and 18. For ease of notation, let

For the sake of simplicity, in Theorem 19 below we have omitted

Let $QB_{k}$ be a rank-$k$ approximation computed by Algorithm 2.2. Then

Remark 3.2 applies to Theorem 20 as well.

Proof of Theorem 20: We only prove Theorem 20 for the Frobenius norm. The case for the 2-norm is completely analogous. As in the proof for Theorem 19, this one involves taking expectations over $\widehat{\Omega}_{2}$ first and $\widehat{\Omega}_{1}$ next. Let $\displaystyle\widehat{\delta}_{k+1}=\sqrt{\sum_{j=k+1}^{n}\sigma^{2}_{j}}$. Fixing $\widehat{\Omega}_{1}$ in Theorem 13 and taking expectation on $\widehat{\Omega}_{2}$ according to Proposition 17, we obtain immediately

For $p\geq 2$, we further take expectation over $\widehat{\Omega}_{1}$ according to Proposition 18. By equation (48),

which is the Frobenius norm upper bound in Theorem 20.

For $p=1$, we again take expectation over $\widehat{\Omega}_{1}$ in equation (53) according to Proposition 18:

which is bounded above by the corresponding expression in Theorem 20.

Now, we turn our attention to the case $p=0$. Taking expectations as before,

Plugging in the expressions for $\alpha$ and $\gamma$ in equation (53),

which is bounded above by the corresponding expression in Theorem 20 since $\widehat{\delta}_{k+1}\leq\sqrt{n-k}\,\,{\sigma_{1}}$. Q.E.D.

3 Large Deviation Bounds

In this section we develop approximation error tail bounds. Theorems 12 and 13 dictate that our main focus will be in developing probabilistic upper bounds on $\left\|\widehat{\Omega}_{2}\right\|_{2}\left\|\widehat{\Omega}_{1}^{{\dagger}}\right\|_{2}$.

with exception probability at most ${\Delta}$.

Proof of Theorem 21: Since $\widehat{\Omega}_{2}$ and $\widehat{\Omega}_{1}$ are independent from each other, we can study how the error depends on the matrix $\widehat{\Omega}_{2}$ when $\widehat{\Omega}_{1}$ is reasonably bounded. To this end, we define an event as follows:

Invoking the conclusion of Lemma 15, we find that

In other words, we have just shown that ${\displaystyle\left\|\widehat{\Omega}_{1}^{{\dagger}}\right\|_{2}\leq t{\cal L}}$ with probability at least $1-t^{-(p+1)}$.

where $X$ has the same dimensions as $\widehat{\Omega}_{2}$. It is straightforward to show that

under event ${\bf E}_{t}$. Also under event ${\bf E}_{t}$ and by Proposition 14, we have

Applying the concentration of measure equation, Theorem 16, conditionally to $\widehat{\Omega}_{2}$ under event ${\bf E}_{t}$,

Use the equation (54) to remove the restriction on $\widehat{\Omega}_{1}$, therefore,

so that $t^{-(p+1)}+e^{-u^{2}/2}=\Delta$. With this choice of $t$ and $u$,

Plugging this bound into the formulas in Theorem 12 and Remark 4.3 proves Theorem 21. Q.E.D.

While the value of oversampling size $p$ does not look so important in the average case error bounds as long as $p\geq 2$, it makes an oversized difference in large deviation bounds. Consider the case $p=2$ with a tiny $\Delta>0$. In this case, ${\displaystyle{\cal C}_{\Delta}}$ may still be quite large, and quite a few extra number of iterations might be necessary to ensure satisfactory convergence with small exception probability.

For $p\leq 1$, the large deviation bound is brutal. For very small values of $k$, such as $1$ in the case of the randomized power method (see Algorithm A.3), it seems unreasonable to require a relatively large value of $p$. On the other hand, a small $p$ value would significantly impact convergence. We will address this conflicting issue of choosing $p$ further in Section 8.

But for any large enough values of $k$ (such as $k=20$ or more, for example,) a reasonable choice would be to choose $p$ so ${\displaystyle\left(\frac{2}{{\Delta}}\right)^{1/(p+1)}}$ is a modest number. We will now choose

This choice gives ${\displaystyle\left(\frac{2}{{\Delta}}\right)^{1/(p+1)}\leq 10}$. For a typical choice of ${\Delta}=10^{-16}$, equation (55) gives $p=16$. For this value of ${\Delta}$, the exception probability is smaller than that of matching DNA fingerprints . Given that the ”random numbers” generated on modern computers are really only pseudo random numbers that may have quite different upper tail distributions than the true Gaussian (see, for example ), and given that only finite precision computations are typically done in practice, it is probably meaningless to require ${\Delta}$ to be much less than $10^{-16}$, the double precision. Additionally, with this choice of $p$, the large deviation bounds are very similar to the average case error bounds, suggesting that the typical behavior is also the worst case behavior, with probability $1-\Delta$.

Our final observation on Theorem 21 is so important that we present it in the form of a Corollary. We will not prove it because it is a direct consequence.

In the notation of Theorem 21, we must have for $j=1,\cdots,k$,

with exception probability at most ${\Delta}$.

This is a surprisingly strong result. We will discuss its implications in terms of rank-revealing factorizations in Section 6 and condition number estimation in Section 7.

Rank-revealing Factorizations

Rank-revealing factorizations were first discussed in Chan . Generally speaking, there are rank-revealing UTV factorizations , QR factorizations , and LU factorizations . While there is no uniform definition of the rank-revealing factorization, a comparison of different forms of rank-revealing factorizations has appeared in Foster and Liu . For the discussions in this section, we make the following definition, which is loosely consistent with those in .

Given $m\times n$ matrices $A$ and $B$ and integer $k<\min(n,m)$, we call $B$ a rank-revealing rank-$k$ approximation to $A$ if ${\bf rank}(B)\leq k$ and if there exist polynomials $c_{1}(m,n)$, and $c_{2}(m,n)$ such that

A rank-revealing rank-$k$ approximation differs from an ordinary rank-$k$ approximation in the extra condition (57), which requires some accuracy in all $k$ leading singular values. Therefore a rank-revealing rank-$k$ approximation is likely a stronger approximation than a simple low rank approximation. To see why (57) is so important, we consider for an example the case where the leading $k+1$ singular values of $A$ are identical: $\sigma_{1}(A)=\cdots=\sigma_{k+1}(A)$. This includes the $n\times n$ identity matrix as a special case. Now choose $\theta=1$ in equation (4). It follows that $B=0$ is an optimal rank-$k$ approximation to $A$, which is likely unacceptable to most users. On the other hand, $B=0$ obviously does not satisfy condition (57) for any polynomial $c_{2}(m,n)$, and therefore is not a rank-revealing rank-$k$ approximation to $A$. Similarly, any orthogonal matrix $Q$ would satisfy the bound in Theorem 2 for such an $A$ matrix, and only the matrix $Q$ from Algorithm 2.2 would satisfy Theorem 21.

By definition, Algorithm 2.2 produces a rank-revealing rank-$k$ approximation with probability at least $1-\Delta$. In this section, we compare this approximation with the strong RRQR factorization developed in Gu and Eisenstat .

(Gu and Eisenstat ) Let $A$ be an $m\times n$ matrix and let $1\leq k\leq\min(m,n)$. For any given parameter $f>1$, there exists a permutation $\Pi$ such that

where for any $1\leq i\leq k$ and $1\leq j\leq n-k$,

Let {\displaystyle\widehat{A}_{k}=Q\left(\begin{array}[]{cc}R_{11}&R_{12}\\ &0\end{array}\right)\Pi^{T}}. Then $\widehat{A}_{k}$ is a rank-$k$ matrix. It follows from equation (59) that

These properties are compatible with the inequalities in Theorem 21. The strong RRQR factorization in Theorem 23 also includes a permutation $\Pi$ that selects $k$ linearly independent columns of $A$ such that $\left\|R_{11}^{-1}R_{12}\right\|_{2}\leq f$. Such information could be useful in some applications .

But the matrix $QB_{k}$, being a two-sided orthogonal approximation, does not contain any information about such permutation. On the other hand, it is likely to be cheaper to compute due to the matrix-matrix product operations involved, and for rapidly decaying singular values or by potentially increasing the value of $q$, it could make a much better approximation than $\widehat{A}_{k}$.

Condition Number Estimation

For any given square non-singular matrix $A$, define

as its condition number. Here $\|\cdot\|$ is any matrix norm, such as the matrix $1$-norm, $2$-norm, $\infty$-norm, Frobenius norm, or $\max$-norm. Condition numbers are of central importance in solving many matrix computation problems, such as linear equations, least squares problems, eigenvalue/eigenvector problems, and sparse matrix problems. For a detailed discussion of condition number estimation, see the survey paper by Higham and the references therein. More recent work includes Laub and Xia .

A typical condition estimator uses a matrix norm estimator to estimate $\|A\|$ and $\|A^{-1}\|$ separately, and multiply them together to get an estimate for $\kappa(A)$. A typical matrix norm estimator, in turn, only accesses the matrix $A$ through matrix-matrix or matrix-vector multiplications, without the need to directly access entries of $A$. Thus the costs of estimating $\|A\|$ and $\|A^{-1}\|$ are similar if a factorization for $A$ is available. The goal in matrix norm estimation is to compute a reliable estimate of $\|A\|$ up to a factor that does not grow too fast with the dimension of $A$, perhaps without direct access to entries of $A$, at a cost that is considerably less than that of matrix factorization or inversion, something that is believed to be impossible (see Remark 5.9.)

Compute $y=Ax,\quad z=A^{T}{\bf sign}(y)$.

$x=e_{j},\quad\mbox{where}\quad j={\bf argmax}_{k}|z_{k}|.$

The $e_{j}$ is the $j$-th unit vector. While it could occasionally take much longer, Hager’s method typically takes very few (less than $5$) iterations to converge to a local maximum that is within a reasonable factor (like $10$ or less) of $\|A\|_{1}$. As Algorithm 2.2 already computes a reliable estimate for $\|A\|_{2}$, it is straightforward to combine Algorithms 2.2 and 7.1 to obtain a reliable estimate for $\|A\|_{1}$, which satisfies $\|A\|_{1}\geq\|A\|_{2}/\sqrt{n}$.

Compute rank-$1$ approximation $QB_{1}$ to $A$ using Algorithm 2.2

Set $\widehat{u}$ to be the right singular vector of $QB_{1}$.

Run Algorithm 7.1 on $A$ with initial vector $x=\widehat{u}/\|\widehat{u}\|_{1}$.

Since $QB_{1}$ is a rank-1 matrix, $\widehat{u}$ is straightforward to compute. The number of iterations in Algorithm 7.1 can be restricted to as few as $1$ or $2$. This is because Algorithm 7.1 is only used to find a column whose vector $1$-norm provides the estimate for $\|A\|_{1}$, no local maximum to problem (60) is necessary. Corollary 24 directly follows from Corollary 22.

For any $0<\Delta\ll 1$, the output $\gamma$ from Algorithm 7.2 must satisfy

with exception probability at most ${\Delta}$.

Hager’s method has been generalized by Higham to estimate the matrix $p$-norm for any $p\geq 1$ and the mixed matrix norm $\|A\|_{\alpha,\beta}$ for $\alpha\geq 1$ and $\beta\geq 1$. In particular, the $\max$-norm is the special case with $\alpha=\infty$ and $\beta=1$. Algorithm 7.2 can be trivially generalized to those cases as well, by replacing Hager’s method in Algorithm 7.2 with its generalized version, leading to a Corollary 24-like conclusion for reliability. We omit the details.

Kuczyński and Woźniakowski developed probabilistic error bounds for estimating the condition number using the Lanczos algorithm for unit start vectors under the uniform distribution. However, our results appear to be much stronger.

Below, we demonstrate the robustness of Algorithm 7.2 through the following example. Let

where $\alpha>0,\rho>0$ are scalars, $b>0$ is an $n-1$ dimensional vector, and $\widehat{A}$ is an $(n-1)\times(n-1)$ matrix. If we take the initial vector $x$ in Algorithm 7.1 to be the vector of all $1$’s (the default choice in LAPACK), then Algorithm 7.1 will always return $\alpha+\|b\|_{1}$ as the $1$-norm estimate, regardless of $\rho\widehat{A}$.

Numerical Experiments

In this section we perform numerical experiments to shed more light on randomized algorithms. Our main purpose of these experiments is to provide numerical support to our probabilistic analysis and to demonstrate that different applications can lead to different singular value distributions in the matrix and impose different accuracy requirements, and thus demand different levels of computational effort on the randomized algorithms.

In the case of a small $k$, it seems unreasonable to require a potentially large value of $p$ as suggested in equation (55). However, for a truely small value of $p$, going random is still not enough to overcome the potential problem of slow convergence associated with a poor start matrix in Algorithm 2.2, and some additional work maybe needed (see Sections 5.)

This discussion is particularly relevant for $k=1$, which corresponds to the classical power method, Algorithm A.2, and its randomized version, Algorithm A.3, in Appendix A. Any value of $p>0$ seems to be too much work, but $p=0$ does not lead to fast enough convergence.

According to Corollary 22, Algorithm 2.2 can already compute order of magnitude approximations to all the leading singular values with $q=0$. Thus, an obvious improvement of Algorithm 2.2 for small values of $k$ would be to compute $\Omega$ with Algorithm 1.1 and then compute a subspace approximation with Algorithm 2.1. Algorithm 8.1 below is designed for subspace computations where ${\displaystyle k=O\left(\lceil\log_{10}\left(\frac{2}{{\Delta}}\right)\rceil\right)}$ or smaller.

Improved Randomized Subspace Iteration for small $k$

Set $\Omega$ to be approximate right singular vector matrix.

We perform our experiments with $4000\times 4000$ matrices of the form

where $\{X_{i}\}$ are $n$-dimensional Gaussian random variables with mean and standard deviation $1$, and where $\{Y_{j}\}$ are $n$-dimensional Gaussian random variables with mean $\mu$ and standard deviation $1$. We choose different $\mu$ values to control the ratio of the two leading singular values of $A$.

For the case of large $\sigma_{2}/\sigma_{1}$ ratio, Algorithm 8.1 converged to $\|A\|_{2}$ in about $250$ steps, as opposed to about $350$ steps for Algorithm A.3. For the case of a small $\sigma_{2}/\sigma_{1}$ ratio, both algorithms performed equally well. Algorithm 8.1 converged slightly more quickly, but that is offset by the extra work needed to compute the initial $\Omega$.

Figure 1 confirms our analysis. At the cost of the initial step to obtain a good start vector, Algorithm 8.1 can converge significantly faster than Algorithm A.3.

2 low-rank approximation

In this experiment, we consider a $4000\times 4000$ matrix of the form

where $\{X_{i}\}$ are equi-spaced points on the edge of the disc $\|X-\pmatrix{-1\cr-1}\|_{2}=\sqrt{2}$ and $\{Y_{j}\}$ equi-spaced points on the edge of the disc $\|Y-\pmatrix{2\cr 2}\|_{2}=2\sqrt{2}$ (see Figure 2.) We compare the performance of Algorithms 1.1 and 2.2 against that of svds, the matlab version of ARPACK for finding a few selected singular values of large matrices. We choose $k=50$. The results are summarized in Table 1.

Since the singular values of this matrix decay relatively quickly, Algorithm 1.1 seems to out-perform Algorithm 2.2 for any values of $q>0$. Algorithm 1.1 also outperforms svds. As Algorithm 1.1 mostly computes matrix-matrix products whereas each step of svds involves a matrix-vector product, we would expect Algorithm 1.1 to have even better performance than svds on modern serial and parallel architectures. This example demonstrates that for matrices with fast decaying singular values, randomized algorithms can be as competitive as the best methods for computing highly accurate low-rank approximations.

3 Structured Matrix Computations

In this example, we demonstrate the effectiveness of randomized algorithms for low-rank approximation in the context of structured matrix computations. ${\tt G}3_{\tt circuit}$ is a $1585478\times 1585478$ sparse SPD matrix arising from circuit simulations. It is publicly available in the University of Flordia Sparse Matrix Collection . Figure 3 depicts its sparsity pattern in the symmetric minimum degree ordering . A direct factorization of this matrix creates a large amount of fill-in. In particular, the Schur complement of the leading $1582178\times 1582178$ principal submatrix, to be called $A$, is a $3300\times 3300$ dense submatrix. Here we compute hierarchical semiseparable (HSS) preconditioners to $A$ with the techniques in and report the numbers of preconditioned conjugate gradient (PCG) steps to iteratively solve for a linear system of equations $Ax=b$ for a random right hand side $b$. The PCG is a very popular technique for solving large SPD systems of equations . We refer the reader to for details about the HSS matrix structure and its numerical construction, but emphasize that the key and most time-consuming step for computing HSS preconditioners is to approximate various off-diagonal blocks of the matrix $A$ by matrices of rank $k$ or less. We choose convergence tolerance $\delta=10^{-12}$. The conjugate gradient method (CG) without any preconditioning takes $878$ iterations to reduce the residual below this tolerance.

Table 2 summarizes our results. We can see that all choices of $p$ drastically decrease the number of CG iterations. Howver, the additional reduction in the number of CG iterations is typically small for higher values of $p$. Considering the extra cost involved in higher $p$ values in the construction of HSS preconditioners, it seems that higher $p$ values are ineffective for this application. This example suggests that for the purpose of constructing preconditioners in structured matrix computations, a small $p$ value is typically sufficient to develop highly effective preconditioners. This is consistent with the rule of thumb that typically randomized algorithms require very little oversampling and a value of $p$ in between $10$ to $20$ suffices . In fact, the SVD truncation in Algorithm 2.2 and Algorithm 2.1 is unnecessary for this example.

4 Eigenfaces

Eigenfaces is a well studied method of face recognition based on principal component analysis (PCA), popularised by the seminal work of Turk and Pentland . For more recent work and survey, see and the references therein. In this experiment we demonstrate the effects of randomized algorithms on face recognition.

Typical face recognition starts with a data base of training images, which are then processed as follows:

Calculate the mean of the training images.

Subtract the mean from the training images, obtaining the mean-shifted images.

Calculate a truncated SVD of the mean-shifted images.

Project the mean-shifted images into the singular vector space using the retained singular vectors, obtaining feature vectors.

To classify a new face, one does the following calculations:

Subtract the mean from the new image, obtaining the mean-shifted image.

Project the mean-shifted image into the singular vector space, obtaining a new feature vector.

Find the feature vector in the data base that best matches the new feature vector.

Our face data are obtained from the Database of Faces maintained at the AT&T Laboratories Cambridge . All faces are greyscale images with a consistent resolution. There are ten different images of each of 40 distinct subjects. The size of each image is $92\times 112$ pixels, with $256$ grey levels per pixel. We use 200 of these images, 5 from each individual, as training images, and the remaining ones for classification.

In Figure 4, the first row are the original face images; the second row are eigenfaces with a rank-$10$ truncated SVD, and the third row eigenfaces with a rank-$20$.

In addition to the exact truncated SVD, we also perform image training and classification using Algorithm 1.1 with different $p$ values. The results are summarized in Table 3. It is clear that smaller $p$ values give worse results than truncated SVD, but $p=40$ gives results that are very similar to truncated SVD, even though some of the singular values are accurate to only within $1$ to $2$ digits. This example demonstrates that limited accuracy that goes beyond being correct to within a constant factor is sufficient for some applications.

Conclusions and Future Work

We have presented some interesting results on randomized algorithms within the framework of the subspace iteration method for singular value and low-rank matrix approximations. While randomized algorithms have been primarily considered as an efficient tool to compute low-rank approximations, our results further suggest that they actually compute the much stronger rank-revealing factorizations, and can be used to reliably estimate condition numbers. We have also presented numerical experimental results that support our analysis.

Acknowledgments. The author would like to thank Shengguo Li, Michael Mahoney, Vladimir Rokhlin, Mark Tygert, Jianlin Xia and Chao Yang for many helpful discussions on this subject. He would especially like to thank Joel Tropp, whose interesting talk at UC Berkeley in the Spring of 2010 sparked the author’s interest on the subject that eventually led to this work, and Chris Melgaard, with whom he had extensive discussions about the material presented in this work. Finally, the author would like to thank the anonymous referees who go out of their ways to provide numerous helpful suggestions that greatly improved the presentation of this paper, including a shorter proof for Theorem 8.

Appendix. For numerical stability, Algorithm A.1 below is often performed once every few iterations in subspace iteration methods, to balance efficiency and numerical stability (see Saad .)

Compute $Y=A\Omega$, and QR factorize $QR=Y$.

Below is the classical power method for computing the $2$-norm of a given matrix.

Compute $Y=\left(AA^{T}\right)^{q}A\,\Omega$.

Compute an orthogonal column basis $Q$ for $Y$.

In situations where no useful information about the leading right singular vector is available, the vector $\Omega$ in Algorithm A.2 can also be chosen to be random, to enhance convergence, leading to

Draw a random $n\times 1$ vector $\Omega$.

Compute $Y=\left(AA^{T}\right)^{q}A\,\Omega$.

Compute an orthogonal column basis $Q$ for $Y$.

Appendix S1 Introduction

In the interest of reducing the length of the original paper, we have put some of the non-essential material here. This Supplemental Material is organized as follows: In Section S2 we discuss how the decaying rates of the singular values can affect parameter choices in the randomized algorithms; in Section 8 we present additional supporting numerical experimental results; in Section S4 we provide proofs for the two propositions in the original paper; and in Section S5 we list the facts we have used from calculus.

Appendix S2 Further Convergence Considerations

as the key factor that controls singular value convergence.

First we consider the model where the singular values (except for the first few) decay and satisfy the following equation

for some constant $\alpha>0$ and any $s,t>1$. This model is satisfied when the singular values of $A$ decay exponentially or faster. We wish to show that Algorithm 1.1 performs better than Algorithm 2.2 with $q>0$.

according to the singular value decay model (S2.61). On the other hand, for Algorithm 1.1, the ratio is

which is a tighter upper bound. This comparison suggests that in general there is little convergence advantage of Algorithm 2.2 over Algorithm 1.1 when the singular values decay exponentially or faster.

S2.2 Slowly decaying singular value distributions

Below we consider the model where the singular values (except for the first few) decay and satisfy the following equation

S2.3 Adaptive Randomized Algorithms

From the two different singular value distributions discussed above, it is clear that much research is needed to design an efficient algorithm that can automatically choose the right set of parameters for different singular value distributions within the framework of Algorithm 2.2.

In this section, we will limit our scope and present an adaptive version of Algorithm 2.2, with the assumption that the singular values decay slowly. Our goal is to quickly compute rank-$k$ approximations up to the tolerance provided. This algorithm is motivated by similar work in and will be used later on in our numerical experiments in Section S3.

: Adaptive Randomized Subspace Iteration Method

Compute $Y=\left(AA^{T}\right)^{q}A\,\Omega$.

Compute an orthogonal column basis $Q$ for $Y$.

Compute the SVD of $B$ and the rank-$k$ truncated SVD $B_{k}$.

Quit. Sampling size exceeding limit for the given tolerance

Update $Y=[Y\quad\left(AA^{T}\right)^{q}A\,\Omega]$.

Update the orthogonal column basis $Q$ for $Y$.

Update the SVD of $B$ and the rank-$k$ truncated SVD $B_{k}$.

Appendix S3 Numerical Experiments

. In this section we report more numerical experimental results to shed more light on randomized algorithms.

Latent Semantic Indexing (LSI) is a massive data processing application based on low-rank approximations . A data base of terms and documents is processed to generate a term-document matrix, where each column is a document with each non-zero in the column represents the weighted number of matches to a particular term.

Given a set of terms (a query), LSI attempts to find the document that best matches it in some semantical sense. To do so, LSI computes a rank-$k$ truncated SVD of the term-document matrix so that $A\approx U_{k}S_{k}V_{k}^{T}$.

For any query vector $q$, compute the feature vector $d=\left(q^{T}U\right)S_{k}^{-1}$. The document that most matches $q$ is the row of $V_{k}$ that is the most parallel to $d$.

We use the TDT2 text data . The TDT2 corpus consists of data collected during the first half of 1998 and taken from 6 sources, including 2 newswires (APW, NYT), 2 radio programs (VOA, PRI) and 2 television programs (CNN, ABC). It consists of 11201 on-topic documents which are classified into 96 semantic categories. What is available at is a subset of this corpus, with a total of 9,394 documents and over $36000$ terms.

We performed $1000$ random queries with the truncated SVD for different values of $k$. Then we repeat the same queries with the low-rank approximation computed by Algorithm S2.1 for $q=0,2,4$ and a decreasing set of $\tau$ values. For each $q$ and $\tau$, Algorithm S2.1 automatically stops once $500$ column samples have been reached in computing the low-rank approximation.

Table 4 clearly indicates that better accuracy in randomized algorithms leads to more agreement with the truncated SVD in terms of query matches. Due to the nature of this experiment, an agreement does not always mean a better match. However, Table 4 does give some indication that better accuracy in the low-rank approximation is probably better for LSI. Since $q=2$ looks significantly better than $q=0$, this example indicates that for LSI, it may be necessary to use Algorithm S2.1 with a small but positive $q$ value for best performance.

Appendix S4 Proofs of Propositions 17 and 18

We begin with the following probability tool.

(Chen and Dongarra ) Let $G$ be an $m\times n$ standard Gaussian random matrix with $m\leq n$, and let $f(x)$ denote the probability density function of $\|G^{\dagger}\|_{2}^{-2}$, then $f(x)$ satisfies:

The following classical result, the law of the unconscious statistician, will be very helpful to our analysis.

Let $g(\cdot)$ be a non-negative continuously differentiable function with $g(0)=0$, and let $G$ be a random matrix, we have

We also need to define the following functions

where $\alpha>0$ and $\delta>0$ are constants to be specified later on. It is easy to see tht $g(0)=0$ and $\widehat{g}(0)=0$, and

Proof of Proposition 17: Define a function $h(G)=\|G\|_{2}$. Then by Proposition 14, we have

$h$ is a Lipschitz function on matrices with Lipschitz constant ${\cal L}=1$ (see Theorem 16):

For equation (35), we can rewrite, by way of function $g(x)$ in (S4.63) and Proposition S4.26,

where in the last equation we have used the fact that ${\displaystyle\int^{\infty}_{0}e^{-u^{2}/2}du=\sqrt{\pi/2}}$ and that ${\displaystyle\int^{\infty}_{0}ue^{-u^{2}/2}du=1}$.

Comparing equations (35) and (S4.65), it is clear that we need to seek a ${\cal C}>0$ so that

for all values of $\alpha>0$. This is equivalent to

For ${\cal C}>{\cal E}$, the right hand side reaches its maximum as $\alpha$ approaches $\infty$. Hence it suffices to choose ${\cal C}$ such that

For $m\geq 1$ and $n\geq 1$, we have ${\cal E}\geq 5$. The last equation for ${\cal C}$ is easily satisfied when we choose ${\cal C}={\cal E}+4=\sqrt{m}+\sqrt{n}+7$.

We will now take a similar approach to prove equation (36). We rewrite, by way of function $\widehat{g}(x)$ in (S4.63),

Comparing equations (36) and (S4.68), we now must seek a ${\cal C}>0$ so that

for all values of $\alpha>0$. Equivalently,

The right hand side approaches the maximum value as $\gamma$ approaches $\infty$. Hence ${\cal C}$ must satisfy

Again the choice ${\cal C}=\sqrt{m}+\sqrt{n}+7={\cal E}+4$ satisfies this equation. Q.E.D.

The Proof for Proposition 18 will follow a similar track. However, due to the complications with $p\leq 1$, we will seek help from Lemma 15 instead of Theorem 16 to shorten the estimation process.

Proof of Proposition 18: As in the proof of Proposition 17, we can write

Following arguments similar to those in the proof of Proposition 17, we have for a constant $C>0$ to be later determined,

Below we will derive lower bounds on (S4.70) for the three difference cases of $p$ in Proposition 18. For $p\geq 2$, equation (S4.70) can be simplified as

for all values of $\alpha>0$. This condition is very similar to equation (S4.66). Arguments similar to those used to solve (S4.66) lead to

Now we consider the case $p=1$. We rewrite equation (S4.70) in light of equation (S5.75) in S5:

To prove Proposition 18, we just need to find a constant ${\cal C}$ so that

where the asymptotic term $\alpha^{2}\log{\displaystyle\frac{2\sqrt{1+\alpha^{2}C^{2}}}{\alpha C}}$ behaves like ${\displaystyle O\left(\alpha^{2}\log{\displaystyle\frac{1}{\alpha}}\right)}$ when $\alpha$ is tiny and like ${\displaystyle O\left(\alpha\right)}$ when $\alpha$ is very large. Equation (S4.71) is equivalent to

All the extra terms involving the $\log$ function have added much complexity to the above expression. We cut it down with equations (S5.79) and (S5.80) in Appendix S5 by replacing all relevant expressions involving $\psi=\alpha C$ by their corresponding calculus upper bounds. This gives

It is now time to prove equation (48). Our approach for $p\geq 2$ is similar. We rewrite, by way of function $\widehat{g}(x)$ in (S4.63),

Similarly, we seek a ${\cal C}>0$ so that

This last equation is very similar to equation (S4.69), with the only difference being the coefficients in the second term on the left hand side. Thus its solution similarly satisfies

The special cases $p=0$ and $p=1$ lead to some involved calculations with Lemma 15. Instead, we will appeal to Lemma S4.25, an upper bound on the probability density function of smallest eigenvalue of the Wishart matrix $GG^{T}$. It is a happy coincidence that this upper bound is reasonably tight for $p\leq 1$. By Lemma S4.25,

The integral in equation (S4.73) can be bounded as

Below we further simplify equation (S4.74). For $p=1$, the integral in (S4.74) becomes, according to equation (S5.77) in S5:

Replacing the integral in equation (S4.74), and plugging the resulting upper bound into equation (S4.73), we obtain the desired equation (48) for $p=1$.

Finally we consider the case $p=0$. The integral in equation (S4.74) can be rewritten as

where we have used the substitution $x=y^{2}$. Applying the inequality

to both factors in the denominator above, and utilizing the identity (S5.78) from S5, we bound the integral from above as

which leads to the desired equation (48) for $p=0$. Q.E.D.

Appendix S5 Facts from Calculus

Here we list the facts we have used from calculus. Their proofs have been left out, since they do not provide any additional insight into our analysis. We start with $4$ definite integrals:

where $\alpha$, $A,B,C,D$ are all positive constants. We will also list the following inequalities for any $\psi>0$:

References