An algorithm for the principal component analysis of large data sets

Introduction

Principal component analysis (PCA) is among the most popular tools in machine learning, statistics, and data analysis more generally. PCA is the basis of many techniques in data mining and information retrieval, including the latent semantic analysis of large databases of text and HTML documents described in . In this paper, we compute PCAs of very large data sets via a randomized version of the block Lanczos method, summarized in Section 3 below. The proofs in and show that this method requires only a couple of iterations to produce nearly optimal accuracy, with overwhelmingly high probability (the probability is independent of the data being analyzed, and is typically $1-10^{-15}$ or greater). The randomized algorithm has many advantages, as shown in and ; the present article adapts the algorithm for use with data sets that are too large to be stored in the random-access memory (RAM) of a typical computer system.

Computing a PCA of a data set amounts to constructing a singular value decomposition (SVD) that accurately approximates the matrix $A$ containing the data being analyzed (possibly after suitably “normalizing” $A$, say by subtracting from each column its mean). That is, if $A$ is $m\times n$, then we must find a positive integer $k<\min(m,n)$ and construct matrices $U$, $\Sigma$, and $V$ such that

with $U$ being an $m\times k$ matrix whose columns are orthonormal, $V$ being an $n\times k$ matrix whose columns are orthonormal, and $\Sigma$ being a diagonal $k\times k$ matrix whose entries are all nonnegative. The algorithm summarized in Section 3 below is most efficient when $k$ is substantially less than $\min(m,n)$; in typical real-world applications, $k\ll\min(m,n)$. Most often, the relevant measure of the quality of the approximation in (1) is the spectral norm of the discrepancy $A-U\,\Sigma\,V^{\top}$; see, for example, Section 3 below. The present article focuses on the spectral norm, though our methods produce similar accuracy in the Frobenius/Hilbert-Schmidt norm (see, for example, ).

The procedure of the present article works to minimize the total number of times that the algorithm has to access each entry of the matrix $A$ being approximated. A related strategy is to minimize the total number of disk seeks and to maximize the dimensions of the approximation that can be constructed with a given amount of RAM; the algorithm in takes this latter approach.

In the present paper, the entries of all matrices are real valued; our techniques extend trivially to matrices whose entries are complex valued. The remainder of the article has the following structure: Section 2 explains the motivation behind the algorithm. Section 3 outlines the algorithm. Section 4 details the implementation for very large matrices. Section 5 quantifies the main factors influencing the running-time of the algorithm. Section 6 illustrates the performance of the algorithm via several numerical examples. Section 7 applies the algorithm to a data set of interest in biochemical imaging. Section 8 draws some conclusions and proposes directions for further research.

Informal description of the algorithm

In this section, we provide a brief, heuristic description. Section 3 below provides more details on the algorithm described intuitively in the present section.

Suppose that $k$, $m$, and $n$ are positive integers with $k<m$ and $k<n$, and $A$ is a real $m\times n$ matrix. We will construct an approximation to $A$ such that

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 diagonal real $k\times k$ matrix whose entries are all nonnegative, $\|A-U\,\Sigma\,V^{\top}\|_{2}$ is the spectral ($l^{2}$-operator) norm of $A-U\,\Sigma\,V^{\top}$, and $\sigma_{k+1}$ is the $(k+1)$st greatest singular value of $A$. To do so, we select nonnegative integers $i$ and $l$ such that $l\geq k$ and $(i+2)k\leq n$ (for most applications, $l=k+2$ and $i\leq 2$ is sufficient; $\|A-U\,\Sigma\,V^{\top}\|_{2}$ will decrease as $i$ and $l$ increase), and then identify an orthonormal basis for “most” of the range of $A$ via the following two steps:

Using a random number generator, form a real $n\times l$ matrix $G$ whose entries are independent and identically distributed Gaussian random variables of zero mean and unit variance, and compute the $m\times((i+1)l)$ matrix

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

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

Intuitively, the columns of $Q$ in (4) constitute an orthonormal basis for most of the range of $A$. Moreover, the somewhat simplified algorithm with $i=0$ is sufficient except when the singular values of $A$ decay slowly; see, for example, .

Notice that $Q$ may have many fewer columns than $A$, that is, $k$ may be substantially less than $n$ (this is the case for most applications of principal component analysis). This is the key to the efficiency of the algorithm.

Having identified a good approximation to the range of $A$, we perform some simple linear algebraic manipulations in order to obtain a good approximation to $A$, via the following four steps:

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

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

The matrices $U$, $\Sigma$, and $V$ obtained via Steps 1–6 above satisfy (2); in fact, they satisfy the more detailed bound (8) described below.

Summary of the algorithm

In this section, we will construct a low-rank (say, rank $k$) approximation $U\,\Sigma\,V^{\top}$ to any given real matrix $A$, such that

with high probability (independent of $A$), where $m$ and $n$ are the dimensions of the given $m\times n$ matrix $A$, $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)$st greatest singular value of $A$, and $C$ is a constant determining the probability of failure (the probability of failure is small when $C=10$, negligible when $C=100$). In (8), $i$ is any nonnegative integer such that $(i+2)k\leq n$ (for most applications, $i=1$ or $i=2$ is sufficient; the algorithm becomes less efficient as $i$ increases), and $\|A-U\,\Sigma\,V^{\top}\|_{2}$ is the spectral ($l^{2}$-operator) norm of $A-U\,\Sigma\,V^{\top}$, that is,

To simplify the presentation, we will assume that $n\leq m$ (if $n>m$, then the user can apply the algorithm to $A^{\top}$). In this section, we summarize the algorithm; see and for an in-depth discussion, including proofs of more detailed variants of (8).

The minimal value of the spectral norm $\|A-B\|_{2}$, minimized over all rank-$k$ matrices $B$, is $\sigma_{k+1}$ (see, for example, Theorem 2.5.3 in ). Hence, (8) guarantees that the algorithm summarized below produces approximations of nearly optimal accuracy.

To construct a rank-$k$ approximation to $A$, we could apply $A$ to about $k$ random vectors, in order to identify the part of its range corresponding to the larger singular values. To help suppress the smaller singular values, we apply $A\,(A^{\top}\,A)^{i}$, too. Once we have identified “most” of the range of $A$, we perform some linear-algebraic manipulations in order to recover an approximation satisfying (8).

A numerically stable realization of the scheme outlined in the preceding paragraph is the following. We choose an integer $l\geq k$ such that $(i+1)l\leq n-k$ (it is generally sufficient to choose $l=k+2$; increasing $l$ can improve the accuracy marginally, but increases computational costs), and make the following six steps:

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

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

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

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

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

In the present paper, we assume that the user specifies the rank $k$ of the approximation $U\,\Sigma\,V^{\top}$ being constructed. See for techniques for determining the rank $k$ adaptively, such that the accuracy $\|A-U\,\Sigma\,V^{\top}\|_{2}$ satisfying (8) also meets a user-specified threshold.

Variants of the fast Fourier transform (FFT) permit additional accelerations; see , , and . However, these accelerations have negligible effect on the algorithm running out-of-core. For out-of-core computations, the simpler techniques of the present paper are preferable.

The algorithm described in the present section can underflow or overflow when the range of the floating-point exponent is inadequate for representing simultaneously both the spectral norm $\|A\|_{2}$ and its $(2i+1)$st power $(\|A\|_{2})^{2i+1}$. A convenient alternative is the algorithm described in ; another solution is to process $A/\|A\|_{2}$ rather than $A$.

Out-of-core computations

With suitably large matrices, some steps in Section 3 above require either storage on disk, or on-the-fly computations obviating the need for storing all the entries of the $m\times n$ matrix $A$ being approximated. Conveniently, Steps 2, 4, 5, and 6 involve only matrices having $\mathcal{O}((i+1)\,l\,(m+n))$ entries; we perform these steps using only storage in random-access memory (RAM). However, Steps 1 and 3 involve $A$, which has $mn$ entries; we perform Steps 1 and 3 differently depending on how $A$ is provided, as detailed below in Subsections 4.1 and 4.2.

If $A$ does not fit in memory, but we have access to a computational routine that can evaluate each entry (or row or column) of $A$ individually, then obviously we can perform Steps 1 and 3 using only storage in RAM. Every time we evaluate an entry (or row or column) of $A$ in order to compute part of a matrix product involving $A$ or $A^{\top}$, we immediately perform all computations associated with this particular entry (or row or column) that contribute to the matrix product.

2 Computations with storage on disk

If $A$ does not fit in memory, but is provided as a file on disk, then Steps 1 and 3 require access to the disk. We assume for definiteness that $A$ is provided in row-major format on disk (if $A$ is provided in column-major format, then we apply the algorithm to $A^{\top}$ instead). To construct the matrix product in (11), we retrieve as many rows of $A$ from disk as will fit in memory, form their inner products with the appropriate columns of $G$, store the results in $H^{(0)}$, and then repeat with the remaining rows of $A$. To construct the matrix product in (17), we initialize all entries of $T$ to zeros, retrieve as many rows of $A$ from disk as will fit in memory, add to $T$ the transposes of these rows, weighted by the appropriate entries of $Q$, and then repeat with the remaining rows of $A$. We construct the matrix product in (12) similarly, forming $F=A^{\top}\,H^{(0)}$ first, and $H^{(1)}=A\,F$ second. Constructing the matrix products in (13)–(14) is analogous.

Computational costs

In this section, we tabulate the computational costs of the algorithm described in Section 3, for the particular out-of-core implementations described in Subsections 4.1 and 4.2. We will be using the notation from Section 3, including the integers $i$, $k$, $l$, $m$, and $n$, and the $m\times n$ matrix $A$.

For most applications, $i\leq 2$ suffices. In contrast, the classical Lanczos algorithm generally requires many iterations in order to yield adequate accuracy, making the computational costs of the classical algorithm prohibitive for out-of-core (or parallel) computations (see, for example, Chapter 9 in ).

We denote by $C_{A}$ the number of floating-point operations (flops) required to evaluate all nonzero entries in $A$. We denote by $N_{A}$ the number of nonzero entries in $A$. With on-the-fly evaluation of the entries of $A$, the six steps of the algorithm described in Section 3 have the following costs:

Forming $H^{(0)}$ in (11) costs $C_{A}+\mathcal{O}(l\,N_{A})$ flops. Forming any of the matrix products in (12)–(14) costs $2C_{A}+\mathcal{O}(l\,N_{A})$ flops. Forming $H$ in (15) costs $\mathcal{O}(ilm)$ flops. All together, Step 1 costs $(2i+1)\,C_{A}+\mathcal{O}(il(m+N_{A}))$ flops.

Forming $Q$ in (16) costs $\mathcal{O}(i^{2}l^{2}m)$ flops.

Forming $T$ in (17) costs $C_{A}+\mathcal{O}(il\,N_{A})$ flops.

Forming the SVD of $T$ in (18) costs $\mathcal{O}(i^{2}l^{2}n)$ flops.

Forming $U$, $\Sigma$, and $V$ in Step 6 costs $\mathcal{O}(k(m+n))$ flops.

Summing up the costs for the six steps above, and using the fact that $k\leq l\leq n\leq m$, we see that the full algorithm requires

flops, where $C_{A}$ is the number of flops required to evaluate all nonzero entries in $A$, and $N_{A}$ is the number of nonzero entries in $A$. In practice, we choose $l\approx k$ (usually a good choice is $l=k+2$).

2 Costs with storage on disk

We denote by $j$ the number of floating-point words of random-access memory (RAM) available to the algorithm. With $A$ stored on disk, the six steps of the algorithm described in Section 3 have the following costs (assuming for convenience that $j>2\,(i+1)\,l\,(m+n)$):

Forming $H^{(0)}$ in (11) requires at most $\mathcal{O}(lmn)$ floating-point operations (flops), $\mathcal{O}(mn/j)$ disk accesses/seeks, and a total data transfer of $\mathcal{O}(mn)$ floating-point words. Forming any of the matrix products in (12)–(14) also requires $\mathcal{O}(lmn)$ flops, $\mathcal{O}(mn/j)$ disk accesses/seeks, and a total data transfer of $\mathcal{O}(mn)$ floating-point words. Forming $H$ in (15) costs $\mathcal{O}(ilm)$ flops. All together, Step 1 requires $\mathcal{O}(ilmn)$ flops, $\mathcal{O}(imn/j)$ disk accesses/seeks, and a total data transfer of $\mathcal{O}(imn)$ floating-point words.

Forming $Q$ in (16) costs $\mathcal{O}(i^{2}l^{2}m)$ flops.

Forming $T$ in (17) requires $\mathcal{O}(ilmn)$ floating-point operations, $\mathcal{O}(mn/j)$ disk accesses/seeks, and a total data transfer of $\mathcal{O}(mn)$ floating-point words.

Forming the SVD of $T$ in (18) costs $\mathcal{O}(i^{2}l^{2}n)$ flops.

Forming $U$, $\Sigma$, and $V$ in Step 6 costs $\mathcal{O}(k(m+n))$ flops.

In practice, we choose $l\approx k$ (usually a good choice is $l=k+2$). Summing up the costs for the six steps above, and using the fact that $k\leq l\leq n\leq m$, we see that the full algorithm requires

disk accesses/seeks (where $j$ is the number of floating-point words of RAM available to the algorithm), and a total data transfer of

floating-point words (more specifically, $C_{\rm words}\approx 2(i+1)mn$).

Numerical examples

In this section, we describe the results of several numerical tests of the algorithm of the present paper.

We set $l=k+2$ for all examples, setting $i=3$ for the first two examples, and $i=1$ for the last two, where $i$, $k$, and $l$ are the parameters from Section 3 above. We ran all examples on a laptop with 1.5 GB of random-access memory (RAM), connected to an external hard drive via USB 2.0. The processor was a single-core 32-bit 2-GHz Intel Pentium M, with 2 MB of L2 cache. We ran all examples in Matlab 7.4.0, storing floating-point numbers in RAM using IEEE standard double-precision variables (requiring 8 bytes per real number), and on disk using IEEE standard single-precision variables (requiring 4 bytes per real number).

All our numerical experiments indicate that the quality and distribution of the pseudorandom numbers have little effect on the accuracy of the algorithm of the present paper. We used Matlab’s built-in pseudorandom number generator for all results reported below.

In this subsection, we illustrate the performance of the algorithm with the principal component analysis of three examples, including a computational simulation.

For the first example, we apply the algorithm to the $m\times n$ matrix

where $E$ and $F$ are $m\times m$ and $n\times n$ unitary discrete cosine transforms of the second type (DCT-II), and $S$ is an $m\times n$ matrix whose entries are zero off the main diagonal, with

Clearly, $S_{1,1}$, $S_{2,2}$, …, $S_{n-1,n-1}$, $S_{n,n}$ are the singular values of $A$.

For the second example, we apply the algorithm to the $m\times n$ matrix

where $E$ and $F$ are $m\times m$ and $n\times n$ unitary discrete cosine transforms of the second type (DCT-II), and $S$ is an $m\times n$ matrix whose entries are zero off the main diagonal, with

Clearly, $S_{1,1}$, $S_{2,2}$, …, $S_{n-1,n-1}$, $S_{n,n}$ are the singular values of $A$.

Table 1a summarizes results of applying the algorithm to the first example, storing on disk the matrix being approximated. Table 1b summarizes results of applying the algorithm to the first example, generating on-the-fly the columns of the matrix being approximated.

Table 2a summarizes results of applying the algorithm to the second example, storing on disk the matrix being approximated. Table 2b summarizes results of applying the algorithm to the second example, generating on-the-fly the columns of the matrix being approximated.

The following list describes the headings of the tables:

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

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

$k$ is the parameter from Section 3 above; $k$ is the rank of the approximation being constructed.

$t_{\rm gen}$ is the time in seconds required to generate and store on disk the matrix $A$ being approximated.

$t_{\rm PCA}$ is the time in seconds required to compute the rank-$k$ approximation (the PCA) provided by the algorithm of the present paper.

$\varepsilon_{0}$ is the spectral norm of the difference between the matrix $A$ being approximated and its best rank-$k$ approximation.

$\varepsilon$ is an estimate of the spectral norm of the difference between the matrix $A$ being approximated and the rank-$k$ approximation produced by the algorithm of the present paper. The estimate $\varepsilon$ of the error is accurate to within a factor of two with extraordinarily high probability; the expected accuracy of the estimate $\varepsilon$ of the error is about 10%, relative to the best possible error $\varepsilon_{0}$ (see ). The appendix below details the construction of the estimate $\varepsilon$ of the spectral norm of $D=A-U\Sigma V^{\top}$, where $A$ is the matrix being approximated, and $U\Sigma V^{\top}$ is the rank-$k$ approximation produced by the algorithm of the present paper.

For the third example, we apply the algorithm with $k=3$ to an $m\times 1000$ matrix whose rows are independent and identically distributed (i.i.d.) realizations of the random vector

where $w_{1}$, $w_{2}$, and $w_{3}$ are orthonormal $1\times 1000$ vectors, $\delta$ is a $1\times 1000$ vector whose entries are i.i.d. Gaussian random variables of mean zero and standard deviation $0.1$, and $(\alpha,\beta,\gamma)$ is drawn at random from inside an ellipsoid with axes of lengths $a=1.5$, $b=1$, and $c=0.5$, specifically,

with $r$ drawn uniformly at random from $$,$\varphi$drawn uniformly at random from$[0,2\pi]$, and$\theta$drawn uniformly at random from$[0,\pi]$. We obtained$w_{1}$,$w_{2}$, and$w_{3}$by applying the Gram-Schmidt process to three vectors whose entries were i.i.d. centered Gaussian random variables;$w_{1}$,$w_{2}$, and$w_{3}$are exactly the same in every row, whereas the realizations of$\alpha$,$\beta$,$\gamma$, and$\delta$in the various rows are independent. We generated all the random numbers on-the-fly using a high-quality pseudorandom number generator; whenever we had to regenerate exactly the same matrix (as the algorithm requires with$i>0$), we restarted the pseudorandom number generator with the original seed.

Figure 1a plots the inner product (i.e., correlation) of $w_{1}$ in (28) and the (normalized) right singular vector associated with the greatest singular value produced by the algorithm of the present article. Figure 1a also plots the inner product of $w_{2}$ in (28) and the (normalized) right singular vector associated with the second greatest singular value, as well as the inner product of $w_{3}$ and the (normalized) right singular vector associated with the third greatest singular value. Needless to say, the inner products (i.e., correlations) all tend to 1, as $m$ increases — as they should. Figure 1b plots the time required to run the algorithm of the present paper, generating on-the-fly the entries of the matrix being processed. The running-time is roughly proportional to $m$, in accordance with (20).

2 Measured data

In this subsection, we illustrate the performance of the algorithm with the principal component analysis of images of faces.

We apply the algorithm with $k=50$ to the 393,216 $\times$ 102,042 matrix whose columns consist of images from the FERET database of faces described in and , with each image duplicated three times. For each duplicate, we set the values of a random choice of 10% of the pixels to numbers chosen uniformly at random from the integers 0, 1, …, 254, 255; all pixel values are integers from 0, 1, …, 254, 255. Before processing with the algorithm of the present article, we “normalized” the matrix by subtracting from each column its mean, then dividing the resulting column by its Euclidean norm. The algorithm of the present paper required 12.3 hours to process all 150 GB of this data set stored on disk, using the laptop computer with 1.5 GB of RAM described earlier (at the beginning of Section 6).

Figure 2a plots the computed singular values. Figure 2b displays the computed “eigenfaces” (that is, the left singular vectors) corresponding to the five greatest singular values.

While this example does not directly provide a reasonable means for performing face recognition or any other task of image processing, it does indicate that the sheer brute force of linear algebra (that is, computing a low-rank approximation) can be used directly for processing (or preprocessing) a very large data set. When used alone, this kind of brute force is inadequate for face recognition and other tasks of image processing; most tasks of image processing can benefit from more specialized methods (see, for example, , , and ). Nonetheless, the ability to compute principal component analyses of very large data sets could prove helpful, or at least convenient.

An application

In this section, we apply the algorithm of the present paper to a data set of interest in a currently developing imaging modality known as single-particle cryo-electron microscopy. For an overview of the field, see , , and their compilations of references.

The data set consists of 10,000 two-dimensional images of the (three-dimensional) charge density map of the E. coli 50S ribosomal subunit, projected from uniformly random orientations, then added to white Gaussian noise whose magnitude is 32 times larger than the original images’, and finally rotated by 0, 1, 2, …, 358, 359 degrees. The entire data set thus consists of 3,600,000 images, each 129 pixels wide and 129 pixels high; the matrix being processed is 3,600,000 $\times$ $129^{2}$. We set $i=1$, $k=250$, and $l=k+2$, where $i$, $k$, and $l$ are the parameters from Section 3 above. Processing the data set required 5.5 hours on two 2.8 GHz quad-core Intel Xeon x5560 microprocessors with 48 GB of random-access memory.

Figure 3a displays the 250 computed singular values. Figure 3b displays the computed right singular vectors corresponding to the 25 greatest computed singular values. Figure 3c displays several noisy projections, their versions before adding the white Gaussian noise, and their denoised versions. Each denoised image is the projection of the corresponding noisy image on the computed right singular vectors associated with the 150 greatest computed singular values. The denoising is clearly satisfactory.

Conclusion

The present article describes techniques for the principal component analysis of data sets that are too large to be stored in random-access memory (RAM), and illustrates the performance of the methods on data from various sources, including standard test sets, numerical simulations, and physical measurements. Several of our data sets stored on disk were so large that less than a hundredth of any of them could fit in our computer’s RAM; nevertheless, the scheme always succeeded. Theorems, their rigorous proofs, and their numerical validations all demonstrate that the algorithm of the present paper produces nearly optimal spectral-norm accuracy. Moreover, similar results are available for the Frobenius/Hilbert-Schmidt norm. Finally, the core steps of the procedures parallelize easily; with the advent of widespread multicore and distributed processing, exciting opportunities for further development and deployment abound.

Appendix

In this appendix, we describe a method for estimating the spectral norm $\|D\|_{2}$ of a matrix $D$. This procedure is particularly useful for checking whether an algorithm has produced a good approximation to a matrix (for this purpose, we choose $D$ to be the difference between the matrix being approximated and its approximation). The procedure is a version of the classic power method, and so requires the application of $D$ and $D^{\top}$ to vectors, but does not use $D$ in any other way. Though the method is classical, its probabilistic analysis summarized below was introduced fairly recently in and (see also Section 3.4 of ).

Suppose that $m$ and $n$ are positive integers, and $D$ is a real $m\times n$ matrix. We define $\omega^{(1)}$, $\omega^{(2)}$, $\omega^{(3)}$, … to be real $n\times 1$ column vectors with independent and identically distributed entries, each distributed as a Gaussian random variable of zero mean and unit variance. For any positive integers $j$ and $k$, we define

which is the best estimate of the spectral norm of $D$ produced by $j$ steps of the power method, started with $k$ independent random vectors (see, for example, ). Naturally, when computing $p_{j,k}(D)$, we do not form $D^{\top}\,D$ explicitly, but instead apply $D$ and $D^{\top}$ successively to vectors.

Needless to say, $p_{j,k}(D)\leq\|D\|_{2}$ for any positive $j$ and $k$. A somewhat involved analysis shows that the probability that

The probability in (34) tends to 1 very quickly as $j$ increases. Thus, even for fairly small $j$, the estimate $p_{j,k}(D)$ of the value of $\|D\|_{2}$ is accurate to within a factor of two, with very high probability; we used $j=6$ for all numerical examples in this paper. We used the procedure of this appendix to estimate the spectral norm in (8), choosing $D=A-U\,\Sigma\,V^{\top}$, where $A$, $U$, $\Sigma$, and $V$ are the matrices from (8). We set $k$ for $p_{j,k}(D)$ to be equal to the rank of the approximation $U\,\Sigma\,V^{\top}$ being constructed.

For more information, see , , or Section 3.4 of .

Acknowledgements

We would like to thank the mathematics departments of UCLA and Yale, especially for their support during the development of this paper and its methods. Nathan Halko and Per-Gunnar Martinsson were supported in part by NSF grants DMS0748488 and DMS0610097. Yoel Shkolnisky was supported in part by Israel Science Foundation grant 485/10. Mark Tygert was supported in part by an Alfred P. Sloan Research Fellowship. Portions of the research in this paper use the FERET database of facial images collected under the FERET program, sponsored by the DOD Counterdrug Technology Development Program Office.

References