Fast Monte Carlo Algorithms for Matrices II: Computing a Low-Rank Approximation to a Matrix

Abstract

In many applications, the data consist of (or may be naturally formulated as) an m × n matrix A. It is often of interest to find a low-rank approximation to A, i.e., an approximation D to the matrix A of rank not greater than a specified rank k, where k is much smaller than m and n. Methods such as the singular value decomposition (SVD) may be used to find an approximation to A which is the best in a well-defined sense. These methods require memory and time which are superlinear in m and n; for many applications in which the data sets are very large this is prohibitive. Two simple and intuitive algorithms are presented which, when given an m × n matrix A, compute a description of a low-rank approximation D * to A, and which are qualitatively faster than the SVD. Both algorithms have provable bounds for the error matrix A − D *. For any matrix X, let X F and X 2 denote its Frobenius norm and its spectral norm, respectively. In the first algorithm, c columns of A are randomly chosen. If the m × c matrix C consists of those c columns of A (after appropriate rescaling), then it is shown that from C T C approximations to the top singular values and corresponding singular vectors may be computed. From the computed singular vectors a description D * of the matrix A may be computed such that rank(D *) ≤ k and such that A − D * 2 ξ ≤ min D:rank(D)≤k A − D 2 ξ + poly(k, 1/c) A 2 F holds with high probability for both ξ = 2, F. This algorithm may be implemented without storing the matrix A in random access memory (RAM), provided it can make two passes over the matrix stored in external memory and use O(cm + c 2) additional RAM. The second algorithm is similar except that it further approximates the matrix C by randomly sampling r rows of C to form a r × c matrix W. Thus, it has additional error, but it can be implemented in three passes over the matrix using only constant additional RAM. To achieve an additional error (beyond the best rank k approximation) that is at most A 2 F , both algorithms take time which is polynomial in k, 1/ , and log(1/δ), where δ > 0 is a failure probability; the first takes time linear in max(m, n) and the second takes time independent of m and n. Our bounds improve previously published results with respect to the rank parameter k for both the Frobenius and spectral norms. In addition, the proofs for the error bounds use a novel method that makes important use of matrix perturbation theory. The probability distribution over columns of A and the rescaling are crucial features of the algorithms which must be chosen judiciously.

Sampling-based dimension reduction for subspace approximation

“…, a m ), and can be computed in time O(min{mn 2 , m 2 n}) using Singular Value Decomposition (SVD). Some recent work on p = 2 case [1,2,3,4,5,9,12], initiated by a result due to Frieze, Kannan, and Vempala [7], has focused on algorithms for computing a k-dimensional subspace that gives (1 + )-approximation to the optimum in time O(mn·poly(k, 1/ )), i.e., linear in the number of co-ordinates we store. Most of these algorithms, with the exception of [1,12], depend on subroutines that sample poly(k, 1/ ) points from given a 1 , a 2 , .…”

Data Streams: Algorithms and Applications

“…One may have to be innovative in seeking appropriate || · || and f . Other linear algebraic functions are similarly of interest such as estimating eigenvalues, determinants, inverses, matrix multiplication and its applications to maxcut, clustering and other graph problems; see the tomes [77,78,79] for sampling solutions to many of these problems for the case when A is fully stored.…”

Localized bases of eigensubspaces and operator compression

“…In many applications ranging from DNA microarray analysis, facial and object recognition, to web search models, we encounter the following problem (1): Given a large N × N matrix A, one wants to find the best approximation of A by a low-rank matrix D, i.e., min D∈R N×N ;rank ðDÞ≤k ‖A − D‖; [1] where the norm is usually taken as the spectral norm ‖ · ‖ 2 or the Frobenius norm ‖ · ‖ F . The solution to this problem is easily formulated in terms of the singular value decomposition (SVD) of A.…”