A randomized algorithm for principal component analysis

Introduction

Principal component analysis (PCA) is among the most widely used techniques in statistics, data analysis, and data mining. PCA is the basis of many machine learning methods, including the latent semantic analysis of large databases of text and HTML documents described in . Computationally, PCA amounts to the low-rank approximation of a matrix containing the data being analyzed. The present article describes an algorithm for the low-rank approximation of matrices, suitable for PCA. This paper demonstrates both theoretically and via numerical examples that the algorithm efficiently produces low-rank approximations whose accuracies are very close to the best possible.

The canonical construction of the best possible rank-$k$ approximation to a real $m\times n$ matrix $A$ uses the singular value decomposition (SVD) of $A$,

where $U$ is a real unitary $m\times m$ matrix, $V$ is a real unitary $n\times n$ matrix, and $\Sigma$ is a real $m\times n$ matrix whose only nonzero entries appear in nonincreasing order on the diagonal and are nonnegative. The diagonal entries $\sigma_{1}$, $\sigma_{2}$, …, $\sigma_{\min(m,n)-1}$, $\sigma_{\min(m,n)}$ of $\Sigma$ are known as the singular values of $A$. The best rank-$k$ approximation to $A$, with $k<m$ and $k<n$, is

where $\sigma_{k+1}$ is the $(k+1)^{\rm st}$ greatest singular value of $A$. For more information about the SVD, see, for example, Chapter 8 in .

where $\|A-B\|$ is the spectral norm of $A-B$, and $\sigma_{k+1}$ is the $(k+1)^{\rm st}$ greatest singular value of $A$ (see, for example, Chapter 5 in ). Furthermore, the algorithms of require only about $\mathcal{O}(nmk)$ flops to produce a rank-$k$ approximation that (unlike an approximation produced by a pivoted $QR$-decomposition) has been guaranteed to satisfy a bound nearly as strong as (4).

While the accuracy in (4) is sufficient for many applications of low-rank approximation, PCA often involves $m\geq$ 10,000, and a “signal-to-noise ratio” $\sigma_{1}/\sigma_{k+1}\leq 100$, where $\sigma_{1}=\|A\|$ is the greatest singular value of $A$, and $\sigma_{k+1}$ is the $(k+1)^{\rm st}$ greatest. Moreover, the singular values $\leq\sigma_{k+1}$ often arise from noise in the process generating the data in $A$, making the singular values of $A$ decay so slowly that $\sigma_{m}\geq\sigma_{k+1}/10$. When $m\geq$ 10,000, $\sigma_{1}/\sigma_{k+1}\leq 100$, and $\sigma_{m}\geq\sigma_{k+1}/10$, the rank-$k$ approximation $B$ produced by a pivoted $QR$-decomposition algorithm typically satisfies $\|A-B\|\sim\|A\|$ — the “approximation” $B$ is effectively unrelated to the matrix $A$ being approximated! For large matrices whose “signal-to-noise ratio” $\sigma_{1}/\sigma_{k+1}$ is less than 10,000, the $\sqrt{m}$ factor in (4) may be unacceptable. Now, pivoted $QR$-decomposition algorithms are not the only algorithms which can compute a rank-$k$ approximation using $\mathcal{O}(nmk)$ flops. However, other algorithms, such as those of , , , , , , , , , , , , , , , , , , , , , , , , , , , and , also yield accuracies involving factors of at least $\sqrt{m}$ when the singular values $\sigma_{k+1}$, $\sigma_{k+2}$, $\sigma_{k+3}$, … of $A$ decay slowly. (The decay is rather slow if, for example, $\sigma_{k+j}\sim j^{\alpha}\,\sigma_{k+1}$ for $j=1$, $2$, $3$, …, with $-1/2<\alpha\leq 0$. Many of these other algorithms are designed to produce approximations having special properties not treated in the present paper, and their spectral-norm accuracy is good when the singular values decay sufficiently fast. Fairly recent surveys of algorithms for low-rank approximation are available in , , and .)

The algorithm described in the present paper produces a rank-$k$ approximation $B$ to $A$ such that

The algorithm of the present paper is randomized, but succeeds with very high probability; for example, the bound (46) on its accuracy holds with probability greater than $1-10^{-15}$. The algorithm is similar to many recently discussed randomized algorithms for low-rank approximation, but produces approximations of higher accuracy when the singular values $\sigma_{k+1}$, $\sigma_{k+2}$, $\sigma_{k+3}$, … of the matrix being approximated decay slowly; see, for example, or . The algorithm is a variant of that in , and the analysis of the present paper should extend to the algorithm of ; stimulated the authors’ collaboration. The algorithm may be regarded as a generalization of the randomized power methods of and , and in fact we use the latter to ascertain the approximations’ accuracy rapidly and reliably.

The algorithm admits obvious “out-of-core” and parallel implementations (assuming that the user chooses the parameter $i$ in (5) to be reasonably small). As with the algorithms of , , , , , , and , the core steps of the algorithm of the present paper involve the application of the matrix $A$ being approximated and its transpose $A^{\hbox{\scriptsize{\rm T}}}$ to random vectors. The algorithm is more efficient when $A$ and $A^{\hbox{\scriptsize{\rm T}}}$ can be applied rapidly to arbitrary vectors, such as when $A$ is sparse.

Throughout the present paper, we use ${\bf 1}$ to denote an identity matrix. We use ${\bf 0}$ to denote a matrix whose entries are all zeros. For any matrix $A$, we use $\|A\|$ to denote the spectral norm of $A$, that is, $\|A\|$ is the greatest singular value of $A$. Furthermore, the entries of all matrices in the present paper are real valued, though the algorithm and analysis extend trivially to matrices whose entries are complex valued.

The present paper has the following structure: Section 2 collects together various known facts which later sections utilize. Section 3 provides the principal lemmas used in bounding the accuracy of the algorithm in Section 4. Section 4 describes the algorithm of the present paper. Section 5 illustrates the performance of the algorithm via several numerical examples. The appendix, Section 6, proves two lemmas stated earlier in Section 3. We encourage the reader to begin with Sections 4 and 5, referring back to the relevant portions of Sections 2 and 3 as they are referenced.

Preliminaries

In this section, we summarize various facts about matrices and functions. Subsection 2.1 discusses the singular values of arbitrary matrices. Subsection 2.2 discusses the singular values of certain random matrices. Subsection 2.3 observes that a certain function is monotone.

The following trivial technical lemma will be needed in Section 3.

Suppose that $m$ and $n$ are positive integers with $m\geq n$. Suppose further that $A$ is a real $m\times n$ matrix such that the least (that is, the $n^{\it th}$ greatest) singular value $\sigma_{n}$ of $A$ is nonzero.

The following lemma states that the greatest singular value of a matrix $A$ is at least as large as the greatest singular value of any rectangular block of entries in $A$; the lemma is a straightforward consequence of the minimax properties of singular values (see, for example, Section 47 of Chapter 2 in ).

Suppose that $k$, $l$, $m$, and $n$ are positive integers with $k\leq m$ and $l\leq n$. Suppose further that $A$ is a real $m\times n$ matrix, and $B$ is a $k\times l$ rectangular block of entries in $A$.

Then, the greatest singular value of $B$ is at most the greatest singular value of $A$.

The following classical lemma provides an approximation $Q\,S$ to an $n\times l$ matrix $R$ via an $n\times k$ matrix $Q$ whose columns are orthonormal, and a $k\times l$ matrix $S$. As remarked in Observation 4, the proof of this lemma provides a classic algorithm for computing $Q$ and $S$, given $R$. We include the proof since we will be using this algorithm.

Suppose that $k$, $l$, and $n$ are positive integers with $k<l\leq n$, and $R$ is a real $n\times l$ matrix.

Then, there exist a real $n\times k$ matrix $Q$ whose columns are orthonormal, and a real $k\times l$ matrix $S$, such that

where $\rho_{k+1}$ is the $(k+1)^{\it st}$ greatest singular value of $R$.

where $U$ is a real $n\times l$ matrix whose columns are orthonormal, $V$ is a real $l\times l$ matrix whose columns are orthonormal, and $\Sigma$ is a real diagonal $l\times l$ matrix, such that

for $j=1$, $2$, …, $l-1$, $l$, where $\Sigma_{j,j}$ is the entry in row $j$ and column $j$ of $\Sigma$, and $\rho_{j}$ is the $j^{\rm th}$ greatest singular value of $R$. We define $Q$ to be the leftmost $n\times k$ block of $U$, and $P$ to be the rightmost $n\times(l-k)$ block of $U$, so that

We define $S$ to be the uppermost $k\times l$ block of $\Sigma\,V^{\hbox{\scriptsize{\rm T}}}$, and $T$ to be the lowermost $(l-k)\times l$ block of $\Sigma\,V^{\hbox{\scriptsize{\rm T}}}$, so that

Combining (8), (9), (10), (11), and the fact that the columns of $U$ are orthonormal, as are the columns of $V$, yields (7). ∎

In order to compute the matrices $Q$ and $S$ in (7) from the matrix $R$, we can construct (8), and then form $Q$ and $S$ according to (10) and (11). (See, for example, Chapter 8 in for details concerning the computation of the SVD.)

2 Singular values of random matrices

The following lemma provides a highly probable upper bound on the greatest singular value of a square matrix whose entries are independent, identically distributed (i.i.d.) Gaussian random variables of zero mean and unit variance; Formula 8.8 in provides an equivalent formulation of the lemma.

Suppose that $n$ is a positive integer, $G$ is a real $n\times n$ matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and $\gamma$ is a positive real number, such that $\gamma>1$ and

Then, the greatest singular value of $G$ is at most $\sqrt{2n}\,\gamma$ with probability not less than the amount in (12).

Combining Lemmas 2 and 5 yields the following lemma, providing a highly probable upper bound on the greatest singular value of a rectangular matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance.

Suppose that $l$, $m$, and $n$ are positive integers with $n\geq l$ and $n\geq m$. Suppose further that $G$ is a real $l\times m$ matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and $\gamma$ is a positive real number, such that $\gamma>1$ and (12) is nonnegative.

Then, the greatest singular value of $G$ is at most $\sqrt{2n}\,\gamma$ with probability not less than the amount in (12).

The following lemma provides a highly probable lower bound on the least singular value of a rectangular matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance; Formula 2.5 in and the proof of Lemma 4.1 in together provide an equivalent formulation of Lemma 7.

Suppose that $j$ and $l$ are positive integers with $j\leq l$. Suppose further that $G$ is a real $l\times j$ matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and $\beta$ is a positive real number, such that

Then, the least (that is, the $j^{\it th}$ greatest) singular value of $G$ is at least $1/(\sqrt{l}\;\beta)$ with probability not less than the amount in (13).

3 A monotone function

The following technical lemma will be needed in Section 4.

Suppose that $\alpha$ is a nonnegative real number, and $f$ is the function defined on $(0,\infty)$ via the formula

Then, $f$ decreases monotonically for $x>\alpha$.

for any positive real number $x$. The right-hand side of (15) is negative when $x>\alpha$. ∎

Mathematical apparatus

In this section, we provide lemmas to be used in Section 4 in bounding the accuracy of the algorithm of the present paper.

The following lemma, proven in the appendix (Section 6), shows that the product $A\,Q\,Q^{\hbox{\scriptsize{\rm T}}}$ of matrices $A$, $Q$, and $Q^{\hbox{\scriptsize{\rm T}}}$ is a good approximation to a matrix $A$, provided that there exist matrices $G$ and $S$ such that

$Q\,S$ is a good approximation to $(G\,(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}\,A)^{\hbox{\scriptsize{\rm T}}}$, and

there exists a matrix $F$ such that $\|F\|$ is not too large, and $F\,G\,(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}\,A$ is a good approximation to $A$.

Suppose that $i$, $k$, $l$, $m$, and $n$ are positive integers with $k\leq l\leq m\leq n$. Suppose further that $A$ is a real $m\times n$ matrix, $Q$ is a real $n\times k$ matrix whose columns are orthonormal, $S$ is a real $k\times l$ matrix, $F$ is a real $m\times l$ matrix, and $G$ is a real $l\times m$ matrix.

The following lemma, proven in the appendix (Section 6), states that, for any positive integer $i$, matrix $A$, and matrix $G$ whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, with very high probability there exists a matrix $F$ with a reasonably small norm, such that $F\,G\,(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}\,A$ is a good approximation to $A$. This lemma is similar to Lemma 19 of .

Suppose that $i$, $j$, $k$, $l$, $m$, and $n$ are positive integers with $j<k<l<m\leq n$. Suppose further that $A$ is a real $m\times n$ matrix, $G$ is a real $l\times m$ matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and $\beta$ and $\gamma$ are positive real numbers, such that the $j^{\it th}$ greatest singular value $\sigma_{j}$ of $A$ is positive, $\gamma>1$, and

Then, there exists a real $m\times l$ matrix $F$ such that

with probability not less than $\Phi$ defined in (17), where $\sigma_{j}$ is the $j^{\it th}$ greatest singular value of $A$, $\sigma_{j+1}$ is the $(j+1)^{\it st}$ greatest singular value of $A$, and $\sigma_{k+1}$ is the $(k+1)^{\it st}$ greatest singular value of $A$.

Given a matrix $A$, and a matrix $G$ whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, the following lemma provides a highly probable upper bound on the singular values of the product $G\,A$ in terms of the singular values of $A$. This lemma is reproduced from , where it appears as Lemma 20.

Suppose that $j$, $k$, $l$, $m$, and $n$ are positive integers with $k<l$, such that $k+j<m$ and $k+j<n$. Suppose further that $A$ is a real $m\times n$ matrix, $G$ is a real $l\times m$ matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and $\gamma$ is a positive real number, such that $\gamma>1$ and

with probability not less than $\Xi$ defined in (20), where $\rho_{k+1}$ is the $(k+1)^{\it st}$ greatest singular value of $G\,A$, $\sigma_{k+1}$ is the $(k+1)^{\it st}$ greatest singular value of $A$, and $\sigma_{k+j+1}$ is the $(k+j+1)^{\it st}$ greatest singular value of $A$.

The following corollary follows immediately from the preceding lemma, by replacing the matrix $A$ with $(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}\,A$, the integer $k$ with $j$, and the integer $j$ with $k-j$.

Suppose $i$, $j$, $k$, $l$, $m$, and $n$ are positive integers with $j<k<l<m\leq n$. Suppose further that $A$ is a real $m\times n$ matrix, $G$ is a real $l\times m$ matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and $\gamma$ is a positive real number, such that $\gamma>1$ and

with probability not less than $\Psi$ defined in (22), where $\rho_{j+1}$ is the $(j+1)^{\it st}$ greatest singular value of $G\,(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}\,A$, $\sigma_{j+1}$ is the $(j+1)^{\it st}$ greatest singular value of $A$, and $\sigma_{k+1}$ is the $(k+1)^{\it st}$ greatest singular value of $A$.

Description of the algorithm

In this section, we describe the algorithm of the present paper, providing details about its accuracy and computational costs. Subsection 4.1 describes the basic algorithm. Subsection 4.2 tabulates the computational costs of the algorithm. Subsection 4.3 describes a complementary algorithm. Subsection 4.4 describes a computationally more expensive variant that is somewhat more accurate and tolerant to roundoff.

Suppose that $i$, $k$, $m$, and $n$ are positive integers with $2k<m\leq n$, and $A$ is a real $m\times n$ matrix. In this subsection, we will construct an approximation to an SVD of $A$ such that

with very high probability, where $U$ is a real $m\times k$ matrix whose columns are orthonormal, $V$ is a real $n\times k$ matrix whose columns are orthonormal, $\Sigma$ is a real diagonal $k\times k$ matrix whose entries are all nonnegative, $\sigma_{k+1}$ is the $(k+1)^{\rm st}$ greatest singular value of $A$, and $C$ is a constant independent of $A$ that depends on the parameters of the algorithm. (Section 5 will give an empirical indication of the size of $C$, and (46) will give one of our best theoretical estimates to date.)

Intuitively, we could apply $A^{\hbox{\scriptsize{\rm T}}}$ to several random vectors, in order to identify the part of its range corresponding to the larger singular values. To enhance the decay of the singular values, we apply $A^{\hbox{\scriptsize{\rm T}}}\,(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}$ instead. Once we have identified most of the range of $A^{\hbox{\scriptsize{\rm T}}}$, we perform several linear-algebraic manipulations in order to recover an approximation to $A$. (It is possible to obtain a similar, somewhat less accurate algorithm by substituting our short, fat matrix $A$ for $A^{\hbox{\scriptsize{\rm T}}}$, and $A^{\hbox{\scriptsize{\rm T}}}$ for $A$.)

More precisely, we choose an integer $l>k$ such that $l\leq m-k$ (for example, $l=k+12$), and make the following five steps:

Using a random number generator, form a real $l\times m$ matrix $G$ whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and compute the $l\times n$ product matrix

Using an SVD, form a real $n\times k$ matrix $Q$ whose columns are orthonormal, such that there exists a real $k\times l$ matrix $S$ for which

where $\rho_{k+1}$ is the $(k+1)^{\rm st}$ greatest singular value of $R$. (See Observation 4 for details concerning the construction of such a matrix $Q$.)

where $U$ is a real $m\times k$ matrix whose columns are orthonormal, $W$ is a real $k\times k$ matrix whose columns are orthonormal, and $\Sigma$ is a real diagonal $k\times k$ matrix whose entries are all nonnegative. (See, for example, Chapter 8 in for details concerning the construction of such an SVD.)

The following theorem states precisely that the matrices $U$, $\Sigma$, and $V$ satisfy (24). See (46) for a more compact (but less general) formulation.

Suppose that $i$, $k$, $l$, $m$, and $n$ are positive integers with $k<l\leq m-k$ and $m\leq n$, and $A$ is a real $m\times n$ matrix. Suppose further that $\beta$ and $\gamma$ are positive real numbers such that $\gamma>1$,

is nonnegative. Suppose in addition that $U$, $\Sigma$, and $V$ are the matrices produced via the five-step algorithm of the present subsection, given above.

with probability not less than $\Pi$, where $\Pi$ is defined in (32), and $\sigma_{k+1}$ is the $(k+1)^{\it st}$ greatest singular value of $A$.

Observing that $U\,\Sigma\,V^{\hbox{\scriptsize{\rm T}}}=A\,Q\,Q^{\hbox{\scriptsize{\rm T}}}$, it is sufficient to prove that

with probability $\Pi$, where $Q$ is the matrix from (26), since combining (34), (27), (28), and (29) yields (33). We now prove (34).

But, it follows from the fact that the columns of $Q$ are orthonormal that

Combining (36), (37), (38), (35), and (31) yields (34), completing the proof for the case when (35) holds.

For the remainder of the proof, we consider the case when

To prove (34), we will use (16) (which is restated and proven in Lemma 6.19 in the appendix), namely,

for any real $m\times l$ matrix $F$, where $G$ is from (25), and $Q$ and $S$ are from (26). We now choose an appropriate matrix $F$.

First, we define $j$ to be the positive integer such that

where $\sigma_{j}$ is the $j^{\rm th}$ greatest singular value of $A$, and $\sigma_{j+1}$ is the $(j+1)^{\rm st}$ greatest (such an integer $j$ exists due to (39) and the supposition of the theorem that $l\leq m-k$). We then use the matrix $F$ from (18) and (19) associated with this integer $j$, so that (as stated in (18) and (19), which are restated and proven in Lemma 6.21 in the appendix)

with probability not less than $\Phi$ defined in (17). Formula (42) bounds the first term in the right-hand side of (40).

To bound the second term in the right-hand side of (40), we observe that $j\leq k$, due to (41) and the supposition of the theorem that $l\leq m-k$. Combining (26), (25), (23), and the fact that $j\leq k$ yields

with probability not less than $\Psi$ defined in (22). Combining (43) and (44) yields

with probability not less than $\Pi$ defined in (32). The combination of Lemma 8, (30), and the fact that $j\leq k$ justifies the use of $k$ (rather than the $j$ used in (17) for $\Phi$) in the last term in the right-hand side of (32).

Combining (40), (42), (45), (41), (31), and the supposition of the theorem that $l\leq m-k$ yields (34), completing the proof. ∎

Choosing $l=k+12$, $\beta=2.57$, and $\gamma=2.43$ in (32) and (33) yields

with probability greater than $1-10^{-15}$, where $\sigma_{k+1}$ is the $(k+1)^{\rm st}$ greatest singular value of $A$. Numerical experiments (some of which are reported in Section 5) indicate that the factor $100l$ in the right-hand side of (46) is much greater than necessary.

Above, we permit $l$ to be any integer greater than $k$. Stronger theoretical bounds on the accuracy are available when $l\geq 2k$. Indeed, via an analysis similar to the proof of Theorem 13 (using in addition the result stated in the abstract of ), it can be shown that the following six-step algorithm with $l\geq 2k$ produces matrices $U$, $\Sigma$, and $V$ satisfying the bound (33) with its right-hand side reduced by a factor of $\sqrt{l}$:

Using a random number generator, form a real $l\times m$ matrix $G$ whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and compute the $l\times n$ product matrix

Using a pivoted $QR$-decomposition algorithm, form a real $n\times l$ matrix $Q$ whose columns are orthonormal, such that there exists a real $l\times l$ matrix $S$ for which

(See, for example, Chapter 5 in for details concerning the construction of such a matrix $Q$.)

2 Computational costs

In this subsection, we tabulate the number of floating-point operations required by the five-step algorithm described in Subsection 4.1 as applied once to a matrix $A$.

The algorithm incurs the following costs in order to compute an approximation to an SVD of $A$:

Forming $R$ in (25) requires applying $A$ to $il$ column vectors, and $A^{\hbox{\scriptsize{\rm T}}}$ to $(i+1)\,l$ column vectors.

Computing $Q$ in (26) costs $\mathcal{O}(l^{2}\,n)$.

Forming $T$ in (27) requires applying $A$ to $k$ column vectors.

Computing the SVD (28) of $T$ costs $\mathcal{O}(k^{2}\,m)$.

Forming $V$ in (29) costs $\mathcal{O}(k^{2}\,n)$.

Summing up the costs in Steps 1–5 above, and using the fact that $k\leq l\leq m\leq n$, we conclude that the algorithm of Subsection 4.1 costs

floating-point operations, where $C_{A}$ is the cost of applying $A$ to a real $n\times 1$ column vector, and $C_{A^{\hbox{\tiny{\rm T}}}}$ is the cost of applying $A^{\hbox{\scriptsize{\rm T}}}$ to a real $m\times 1$ column vector.

We observe that the algorithm only requires applying $A$ to $il+k$ vectors and $A^{\hbox{\scriptsize{\rm T}}}$ to $il+l$ vectors; it does not require explicit access to the individual entries of $A$. This consideration can be important when $A$ and $A^{\hbox{\scriptsize{\rm T}}}$ are available solely in the form of procedures for their applications to arbitrary vectors. Often such procedures for applying $A$ and $A^{\hbox{\scriptsize{\rm T}}}$ cost much less than the standard procedure for applying a dense matrix to a vector.

3 A modified algorithm

In this subsection, we describe a simple modification of the algorithm described in Subsection 4.1. Again, suppose that $i$, $k$, $l$, $m$, and $n$ are positive integers with $k<l\leq m-k$ and $m\leq n$, and $A$ is a real $m\times n$ matrix. Then, the following five-step algorithm constructs an approximation to an SVD of $A^{\hbox{\scriptsize{\rm T}}}$ such that

with very high probability, where $U$ is a real $n\times k$ matrix whose columns are orthonormal, $V$ is a real $m\times k$ matrix whose columns are orthonormal, $\Sigma$ is a real diagonal $k\times k$ matrix whose entries are all nonnegative, $\sigma_{k+1}$ is the $(k+1)^{\rm st}$ greatest singular value of $A$, and $C$ is a constant independent of $A$ that depends on the parameters of the algorithm:

Using a random number generator, form a real $l\times m$ matrix $G$ whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and compute the $l\times m$ product matrix

Using an SVD, form a real $m\times k$ matrix $Q$ whose columns are orthonormal, such that there exists a real $k\times l$ matrix $S$ for which

where $\rho_{k+1}$ is the $(k+1)^{\rm st}$ greatest singular value of $R$. (See Observation 4 for details concerning the construction of such a matrix $Q$.)

where $U$ is a real $n\times k$ matrix whose columns are orthonormal, $W$ is a real $k\times k$ matrix whose columns are orthonormal, and $\Sigma$ is a real diagonal $k\times k$ matrix whose entries are all nonnegative. (See, for example, Chapter 8 in for details concerning the construction of such an SVD.)

Clearly, (53) is similar to (24), as (54) is similar to (25).

The ideas of Remark 4.15 are obviously relevant to the algorithm of the present subsection, too.

4 Blanczos

In this subsection, we describe a modification of the algorithm of Subsection 4.1, enhancing the accuracy at a little extra computational expense. Suppose that $i$, $k$, $l$, $m$, and $n$ are positive integers with $k<l$ and $(i+1)l\leq m-k$, and $A$ is a real $m\times n$ matrix, such that $m\leq n$. Then, the following five-step algorithm constructs an approximation $U\,\Sigma\,V^{\hbox{\scriptsize{\rm T}}}$ to an SVD of $A$:

Using a random number generator, form a real $l\times m$ matrix $G$ whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and compute the $l\times n$ matrices $R^{(0)}$, $R^{(1)}$, …, $R^{(i-1)}$, $R^{(i)}$ defined via the formulae

Using a pivoted $QR$-decomposition algorithm, form a real $n\times((i+1)l)$ matrix $Q$ whose columns are orthonormal, such that there exists a real $((i+1)l)\times((i+1)l)$ matrix $S$ for which

(See, for example, Chapter 5 in for details concerning the construction of such a matrix $Q$.)

Compute the $m\times((i+1)l)$ product matrix

where $U$ is a real $m\times((i+1)l)$ matrix whose columns are orthonormal, $W$ is a real $((i+1)l)\times((i+1)l)$ matrix whose columns are orthonormal, and $\Sigma$ is a real diagonal $((i+1)l)\times((i+1)l)$ matrix whose entries are all nonnegative. (See, for example, Chapter 8 in for details concerning the construction of such an SVD.)

Compute the $n\times((i+1)l)$ product matrix

An analysis similar to the proof of Theorem 13 above shows that the matrices $U$, $\Sigma$, and $V$ produced by the algorithm of the present subsection satisfy the same upper bounds (33) and (46) as the matrices produced by the algorithm of Subsection 4.1. If desired, one may produce a similarly accurate rank-$k$ approximation by arranging $U$, $\Sigma$, and $V$ such that the diagonal entries of $\Sigma$ appear in nonincreasing order, and then discarding all but the leftmost $k$ columns of $U$ and all but the leftmost $k$ columns of $V$, and retaining only the leftmost uppermost $k\times k$ block of $\Sigma$. We will refer to the algorithm of the present subsection as “blanczos,” due to its similarity with the block Lanczos method (see, for example, Subsection 9.2.6 in for a description of the block Lanczos method).

Numerical results

In this section, we illustrate the performance of the algorithm of the present paper via several numerical examples.

We use the algorithm to construct a rank-$k$ approximation, with $k=10$, to the $m\times(2m)$ matrix $A$ defined via its singular value decomposition

where $U^{(A)}$ is an $m\times m$ Hadamard matrix (a unitary matrix whose entries are all $\pm 1/\sqrt{m}$), $V^{(A)}$ is a $(2m)\times(2m)$ Hadamard matrix, and $\Sigma^{(A)}$ is an $m\times(2m)$ matrix whose entries are zero off the main diagonal, and whose diagonal entries are defined in terms of the $(k+1)^{\rm st}$ singular value $\sigma_{k+1}$ via the formulae

for $j=1$, $2$, …, $9$, $10$, where $\lfloor j/2\rfloor$ is the greatest integer less than or equal to $j/2$, and

for $j=11$, $12$, …, $m-1$, $m$. Thus, $\sigma_{1}=1$ and $\sigma_{k}=\sigma_{k+1}$ (recall that $k=10$). We always choose $\sigma_{k+1}<1$, so that $\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{m-1}\geq\sigma_{m}$.

Figure 1 plots the singular values $\sigma_{1}$, $\sigma_{2}$, …, $\sigma_{m-1}$, $\sigma_{m}$ of $A$ with $m=512$ and $\sigma_{k+1}=.001$; these parameters correspond to the first row of numbers in Table 1, the first row of numbers in Table 2, and the first row of numbers in Table 6.

Table 1 reports the results of applying the five-step algorithm of Subsection 4.1 to matrices of various sizes, with $i=1$. Table 2 reports the results of applying the five-step algorithm of Subsection 4.1 to matrices of various sizes, with $i=0$. The algorithms of , , and for low-rank approximation are essentially the same as the algorithm used for Table 2 (with $i=0$).

Table 3 reports the results of applying the five-step algorithms of Subsections 4.1 and 4.3 with varying numbers of iterations $i$. Rows in the table where $i$ is enclosed in parentheses correspond to the algorithm of Subsection 4.3; rows where $i$ is not enclosed in parentheses correspond to the algorithm of Subsection 4.1.

Table 4 reports the results of applying the five-step algorithm of Subsection 4.1 to matrices whose best rank-$k$ approximations have varying accuracies. Table 5 reports the results of applying the blanczos algorithm of Subsection 4.4 to matrices whose best rank-$k$ approximations have varying accuracies.

Table 6 reports the results of calculating pivoted $QR$-decompositions, via plane (Householder) reflections, of matrices of various sizes. We computed the pivoted $QR$-decomposition of the transpose of $A$ defined in (69), rather than of $A$ itself, for reasons of accuracy and efficiency. As pivoted $QR$-decomposition requires dense matrix arithmetic, our 1 GB of random-access memory (RAM) imposed the limit $m\leq 4096$ for Table 6.

The headings of the tables have the following meanings:

$m$ is the number of rows in $A$, the matrix being approximated.

$n$ is the number of columns in $A$, the matrix being approximated.

$i$ is the integer parameter used in the algorithms of Subsections 4.1, 4.3, and 4.4. Rows in the tables where $i$ is enclosed in parentheses correspond to the algorithm of Subsection 4.3; rows where $i$ is not enclosed in parentheses correspond to either the algorithm of Subsection 4.1 or that of Subsection 4.4.

$t$ is the time in seconds required by the algorithm to create an approximation and compute its accuracy $\delta$.

$\sigma_{k+1}$ is the $(k+1)^{\rm st}$ greatest singular value of $A$, the matrix being approximated; $\sigma_{k+1}$ is also the accuracy of the best possible rank-$k$ approximation to $A$.

$\delta$ is the accuracy of the approximation $U\,\Sigma\,V^{\hbox{\scriptsize{\rm T}}}$ (or $(QRP)^{\hbox{\scriptsize{\rm T}}}$, for Table 6) constructed by the algorithm. For Tables 1–5,

where $U$ is an $m\times k$ matrix whose columns are orthonormal, $V$ is an $n\times k$ matrix whose columns are orthonormal, and $\Sigma$ is a diagonal $k\times k$ matrix whose entries are all nonnegative; for Table 6,

where $P$ is an $m\times m$ permutation matrix, $R$ is a $k\times m$ upper-triangular (meaning upper-trapezoidal) matrix, and $Q$ is an $n\times k$ matrix whose columns are orthonormal.

The values for $t$ are the average values over 3 independent randomized trials of the algorithm. The values for $\delta$ are the worst (maximum) values encountered in 3 independent randomized trials of the algorithm. The values for $\delta$ in each trial are those produced by 20 iterations of the power method applied to $A-U\,\Sigma\,V^{\hbox{\scriptsize{\rm T}}}$ (or $A-(QRP)^{\hbox{\scriptsize{\rm T}}}$, for Table 6), started with a vector whose entries are i.i.d. centered Gaussian random variables. The theorems of and guarantee that this power method produces accurate results with overwhelmingly high probability.

We performed all computations using IEEE standard double-precision variables, whose mantissas have approximately one bit of precision less than 16 digits (so that the relative precision of the variables is approximately .2E–15). We ran all computations on one core of a 1.86 GHz Intel Centrino Core Duo microprocessor with 2 MB of L2 cache and 1 GB of RAM. We compiled the Fortran 77 code using the Lahey/Fujitsu Linux Express v6.2 compiler, with the optimization flag --o2 enabled. We implemented a fast Walsh-Hadamard transform to apply rapidly the Hadamard matrices $U^{(A)}$ and $V^{(A)}$ in (69). We used plane (Householder) reflections to compute all pivoted $QR$-decompositions. We used the LAPACK 3.1.1 divide-and-conquer SVD routine dgesdd to compute all full SVDs. For the parameter $l$, we set $l=12$ $(=k+2)$ for all of the examples reported here.

The experiments reported here and our further tests point to the following:

The accuracies in Table 1 are superior to those in Table 2; the algorithm performs much better with $i>0$. (The algorithms of , , and for low-rank approximation are essentially the same as the algorithm used for Tables 1 and 2 when $i=0$.)

The accuracies in Table 1 are superior to the corresponding accuracies in Table 6; the algorithm of the present paper produces higher accuracy than the classical pivoted $QR$-decompositions for matrices whose spectra decay slowly (such as those matrices tested in the present section).

The accuracies in Tables 1–3 appear to be proportional to $m^{1/(4i+2)}\,\sigma_{k+1}$ for the algorithm of Subsection 4.1, and to be proportional to $m^{1/(4i)}\,\sigma_{k+1}$ for the algorithm of Subsection 4.3, in accordance with (24) and (53). The numerical results reported here, as well as our further experiments, indicate that the theoretical bound (33) on the accuracy should remain valid with a greatly reduced constant in the right-hand side, independent of the matrix $A$ being approximated. See item 6 below for a discussion of Tables 4 and 5.

The timings in Tables 1–5 are consistent with (52), as we could (and did) apply the Hadamard matrices $U^{(A)}$ and $V^{(A)}$ in (69) to vectors via fast Walsh-Hadamard transforms at a cost of $\mathcal{O}(m\,\log(m))$ floating-point operations per matrix-vector multiplication.

The quality of the pseudorandom number generator has almost no effect on the accuracy of the algorithm, nor does substituting uniform variates for the normal variates.

The accuracies in Table 5 are superior to those in Table 4, particularly when the $k^{\rm th}$ greatest singular value $\sigma_{k}$ of the matrix $A$ being approximated is very small. Understandably, the algorithm of Subsection 4.1 would seem to break down when $(\sigma_{k})^{2i+1}$ is less than the machine precision, while $\sigma_{k}$ itself is not, unlike the blanczos algorithm of Subsection 4.4. When $(\sigma_{k})^{2i+1}$ is much less than the machine precision, while $\sigma_{k}$ is not, the accuracy of blanczos in the presence of roundoff is similar to that of the algorithm of Subsection 4.1 run with a reduced $i$. When $(\sigma_{k})^{2i+1}$ is much greater than the machine precision, the accuracy of blanczos is similar to that of the algorithm of Subsection 4.1 run with $i$ being the same as in the blanczos algorithm. Since the blanczos algorithm of Subsection 4.4 is so tolerant of roundoff, we suspect that the blanczos algorithm is a better general-purpose black-box tool for the computation of principal component analyses, despite its somewhat higher cost as compared with the algorithms of Subsections 4.1 and 4.3.

A MATLAB{}^{\hbox{\ooalign{\hfil\raise 0.60275pt\hbox{{\tiny R}}\hfil\crcr\text{\mathchar 524\relax}}}} implementation of the blanczos algorithm of Subsection 4.4 is available on the file exchange at http://www.mathworks.com in the package entitled, “Principal Component Analysis.”

Appendix

In this appendix, we restate and prove Lemmas 9 and 10 from Section 3.

The following lemma, stated earlier as Lemma 9 in Section 3, shows that the product $A\,Q\,Q^{\hbox{\scriptsize{\rm T}}}$ of matrices $A$, $Q$, and $Q^{\hbox{\scriptsize{\rm T}}}$ is a good approximation to a matrix $A$, provided that there exist matrices $G$ and $S$ such that

$Q\,S$ is a good approximation to $(G\,(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}\,A)^{\hbox{\scriptsize{\rm T}}}$, and

there exists a matrix $F$ such that $\|F\|$ is not too large, and $F\,G\,(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}\,A$ is a good approximation to $A$.

Suppose that $i$, $k$, $l$, $m$, and $n$ are positive integers with $k\leq l\leq m\leq n$. Suppose further that $A$ is a real $m\times n$ matrix, $Q$ is a real $n\times k$ matrix whose columns are orthonormal, $S$ is a real $k\times l$ matrix, $F$ is a real $m\times l$ matrix, and $G$ is a real $l\times m$ matrix.

The proof is straightforward, but tedious, as follows.

We obtain from the triangle inequality that

First, we provide a bound for $\|A\,Q\,Q^{\hbox{\scriptsize{\rm T}}}-F\,G\,B\,Q\,Q^{\hbox{\scriptsize{\rm T}}}\|$. Clearly,

It follows from the fact that the columns of $Q$ are orthonormal that

Next, we provide a bound for $\|F\,G\,B\,Q\,Q^{\hbox{\scriptsize{\rm T}}}-F\,G\,B\|$. Clearly,

It follows from the triangle inequality that

Also, it follows from the fact that the columns of $Q$ are orthonormal that

Combining (76), (80), (88), and (75) yields (74).

The following lemma, stated earlier as Lemma 10 in Section 3, shows that, for any positive integer $i$, matrix $A$, and matrix $G$ whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, with very high probability there exists a matrix $F$ with a reasonably small norm, such that $F\,G\,(A\,A^{\hbox{\scriptsize{\rm T}}})^{i}\,A$ is a good approximation to $A$. This lemma is similar to Lemma 19 of .

Suppose that $i$, $j$, $k$, $l$, $m$, and $n$ are positive integers with $j<k<l<m\leq n$. Suppose further that $A$ is a real $m\times n$ matrix, $G$ is a real $l\times m$ matrix whose entries are i.i.d. Gaussian random variables of zero mean and unit variance, and $\beta$ and $\gamma$ are positive real numbers, such that the $j^{\it th}$ greatest singular value $\sigma_{j}$ of $A$ is positive, $\gamma>1$, and

Then, there exists a real $m\times l$ matrix $F$ such that

with probability not less than $\Phi$ defined in (89), where $\sigma_{j}$ is the $j^{\it th}$ greatest singular value of $A$, $\sigma_{j+1}$ is the $(j+1)^{\it st}$ greatest singular value of $A$, and $\sigma_{k+1}$ is the $(k+1)^{\it st}$ greatest singular value of $A$.

We prove the existence of a matrix $F$ satisfying (90) and (91) by constructing one.

where $U$ is a real unitary $m\times m$ matrix, $\Sigma$ is a real diagonal $m\times m$ matrix, and $V$ is a real $n\times m$ matrix whose columns are orthonormal, such that

for $p=1$, $2$, …, $m-1$, $m$, where $\Sigma_{p,p}$ is the entry in row $p$ and column $p$ of $\Sigma$, and $\sigma_{p}$ is the $p^{\rm th}$ greatest singular value of $A$.

Next, we define auxiliary matrices $H$, $R$, $\Gamma$, $S$, $T$, $\Theta$, and $P$. We define $H$ to be the leftmost $l\times j$ block of the $l\times m$ matrix $G\,U$, $R$ to be the $l\times(k-j)$ block of $G\,U$ whose first column is the $(k+1)^{\rm st}$ column of $G\,U$, and $\Gamma$ to be the rightmost $l\times(m-k)$ block of $G\,U$, so that

Combining the fact that $U$ is real and unitary, and the fact that the entries of $G$ are i.i.d. Gaussian random variables of zero mean and unit variance, we see that the entries of $H$ are also i.i.d. Gaussian random variables of zero mean and unit variance, as are the entries of $R$, and as are the entries of $\Gamma$. We define $H^{(-1)}$ to be the real $j\times l$ matrix given by the formula

($H^{\hbox{\scriptsize{\rm T}}}\,H$ is invertible with high probability due to Lemma 7). We define $S$ to be the leftmost uppermost $j\times j$ block of $\Sigma$, $T$ to be the $(k-j)\times(k-j)$ block of $\Sigma$ whose leftmost uppermost entry is the entry in the $(j+1)^{\rm st}$ row and $(j+1)^{\rm st}$ column of $\Sigma$, and $\Theta$ to be the rightmost lowermost $(m-k)\times(m-k)$ block of $\Sigma$, so that

We define $P$ to be the real $m\times l$ matrix whose uppermost $j\times l$ block is the product $S^{-2i}\,H^{(-1)}$, whose entries are zero in the $(k-j)\times l$ block whose first row is the $(j+1)^{\rm st}$ row of $P$, and whose entries in the lowermost $(m-k)\times l$ block are zero, so that

Finally, we define $F$ to be the $m\times l$ matrix given by

Combining (95), (6), the fact that the entries of $H$ are i.i.d. Gaussian random variables of zero mean and unit variance, and Lemma 7 yields

Combining (98), (99), (96), (93), the fact that $\Sigma$ is zero off its main diagonal, and the fact that $U$ is unitary yields (91).

We now show that $F$ defined in (98) satisfies (90).

Combining (101)–(108) and the fact that the columns of $U$ are orthonormal, as are the columns of $V$, yields

Combining Lemma 6 and the fact that the entries of $R$ are i.i.d. Gaussian random variables of zero mean and unit variance, as are the entries of $\Gamma$, yields

Combining (109), (99), (110), and (111) yields

with probability not less than $\Phi$ defined in (89). Combining (113), the fact that $\sigma_{j+1}\leq\sigma_{j}$, and the fact that

for any nonnegative real numbers $x$ and $y$ yields (90).

Acknowledgements

We thank Ming Gu for suggesting the combination of the Lanczos method with randomized methods for the low-rank approximation of matrices. We are grateful for many helpful discussions with R. Raphael Coifman and Yoel Shkolnisky. We thank the anonymous referees for their useful suggestions.

References