We propose a Standing Wave Decomposition (SWD) approximation to Gaussian Process regression (GP). GP involves a costly matrix inversion operation, which limits applicability to large data analysis. For an input space that can be approximated by a grid and when correlations among data are short-ranged, the kernel matrix inversion can be replaced by analytic diagonalization using the SWD. We show that this approach applies to uni- and multi-dimensional input data, extends to include longer-range correlations, and the grid can be in a latent space and used as inducing points. Through simulations, we show that our approximate method outperforms existing methods in predictive accuracy per unit time in the regime where data are plentiful. Our SWD-GP is recommended for regression analyses where there is a relatively large amount of data and or there are constraints on computation time.
Previous works have also taken advantage of structured data to increase the computational efficiency of GP. For input data that lie on a grid, @cite_17 @cite_8 sidestep the matrix inversion by solving @math with gradient-based methods. Importantly, the optimization is fast because Toeplitz circulant structure can be exploited for fast matrix multiplication. These fast multiplications take on forms similar to . For input data that do not lie on a grid, @cite_11 (KISS-GP) introduces a latent grid, finds a sparse, approximate representation of @math in terms of @math , and makes inference using the same kind of optimization problem in @cite_17 @cite_8 . The sparsification in KISS-GP and that in our LGSWD-GP share the same flavor in ignoring the diminishing long-range correlations among data points. Recently, @cite_24 further extended this optimization approach to efficiently handle high dimensional data by applying the Nystr " o m approximation for eigen-decomposition and by exploiting properties of Kronecker and Khatri-Rao products for fast matrix multiplication. Overall, our method and these methods all take advantage of the grid structure to form fast matrix multiplication; however, our method differs in that it overcomes the matrix inversion bottleneck via analytic diagonalizaion rather than optimization.
@cite_24 @cite_8 @cite_11 @cite_17
1252747413 2061781302 2075490698 2035233604 2963396025 2089135543 1990568697 1521197246 2029301183 2124007994 2126407450 2109449402 1712535590 2952294252 2086953401 2045748842 2020461703 2141719776 1522470213 1998269045
This paper presents a fast direct solver for 3D discretized linear systems using the supernodal multifrontal method together with low-rank approximations. For linear systems arising from certain partial differential equations (PDEs) such as elliptic equations, during the Gaussian elimination of the matrices with Nested Dissection ordering, the fill-in of L and U factors loses its sparsity and contains dense blocks with low-rank property. Off-diagonal blocks can be efficiently approximated with low-rank matrices; diagonal blocks approximated with semiseparable structures called hierarchically semiseparable (HSS) representations. Matrix operations in the multifrontal method are performed in low-rank arithmetic. We present efficient way to organize the HSS structured operations along the elimination. To compress dense blocks into low-rank or HSS structures, we use effective cross approximation (CA) approach. We also use idea of adaptive balancing between robust arithmetic for computing the small dense blocks and low-rank matrix operations for handling with compressed ones while performing the Gaussian elimination. This new proposed solver can be essentially parallelized both on architecture with shared and distributed memory and can be used as effective preconditioner. To check efficient of our solver we compare it with Intel MKL PARDISO - the high performance direct solver. Memory and performance tests demonstrate up to 3 times performance and memory gain for the 3D problems with more than (10^6 ) unknowns. Therefore, proposed multifrontal HSS solver can solve large problems, which cannot be resolved by direct solvers because of large memory consumptions. The goal of this paper is to present an efficient implementation of an explicit matrix inversion of general square matrices on multicore computer architecture. The inversion procedure is split into four steps: 1) computing the LU factorization, 2) inverting the upper triangular U factor, 3) solving a linear system, whose solution yields inverse of the original matrix and 4) applying backward column pivoting on the inverted matrix. Using a tile data layout, which represents the matrix in the system memory with an optimized cache-aware format, the computation of the four steps is decomposed into computational tasks. A directed acyclic graph is generated on the fly which represents the program data flow. Its nodes represent tasks and edges the data dependencies between them. Previous implementations of matrix inversions, available in the state-of-the-art numerical libraries, are suffer from unnecessary synchronization points, which are non-existent in our implementation in order to fully exploit the parallelism of the underlying hardware. Our algorithmic approach allows to remove these bottlenecks and to execute the tasks with loose synchronization. A runtime environment system called QUARK is necessary to dynamically schedule our numerical kernels on the available processing units. The reported results from our LU-based matrix inversion implementation significantly outperform the state-of-the-art numerical libraries such as LAPACK (5x), MKL (5x) and ScaLAPACK (2.5x) on a contemporary AMD platform with four sockets and the total of 48 cores for a matrix of size 24000. A power consumption analysis shows that our high performance implementation is also energy efficient and substantially consumes less power than its competitors. In this paper we develop a fast direct solver for large discretized linear systems using the supernodal multifrontal method together with low-rank approximations. For linear systems arising from certain partial differential equations such as elliptic equations, during the Gaussian elimination of the matrices with proper ordering, the fill-in has a low-rank property: all off-diagonal blocks have small numerical ranks with proper definition of off-diagonal blocks. Matrices with this low-rank property can be efficiently approximated with semiseparable structures called hierarchically semiseparable (HSS) representations. We reveal the above low-rank property by ordering the variables with nested dissection and eliminating them with the multifrontal method. All matrix operations in the multifrontal method are performed in HSS forms. We present efficient ways to organize the HSS structured operations along the elimination. Some fast HSS matrix operations using tree structures are proposed. This new structured multifrontal method has nearly linear complexity and a linear storage requirement. Thus, we call it a superfast multifrontal method. It is especially suitable for large sparse problems and also has natural adaptability to parallel computations and great potential to provide effective preconditioners. Numerical results demonstrate the efficiency. In this paper, we study decomposition methods based on separable approximations for minimizing the augmented Lagrangian. In particular, we study and compare the diagonal quadratic approximation method DQAM of Mulvey and Ruszczynski [A diagonal quadratic approximation method for large scale linear programs, Oper. Res. Lett. 12 1992, pp. 205–215] and the parallel coordinate descent method PCDM of Richtarik and Takac [Parallel coordinate descent methods for big data optimization. Technical report, November 2012. arXiv:1212.0873]. We show that the two methods are equivalent for feasibility problems up to the selection of a step-size parameter. Furthermore, we prove an improved complexity bound for PCDM under strong convexity, and show that this bound is at least 8L′ Lω−12 times better than the best known bound for DQAM, where ω is the degree of partial separability and L’ and L are the maximum and average of the block Lipschitz constants of the gradient of the quadratic penalty appearing in the augmented Lagrangian. We introduce an online tensor decomposition based approach for two latent variable modeling problems namely, (1) community detection, in which we learn the latent communities that the social actors in social networks belong to, and (2) topic modeling, in which we infer hidden topics of text articles. We consider decomposition of moment tensors using stochastic gradient descent. We conduct optimization of multilinear operations in SGD and avoid directly forming the tensors, to save computational and storage costs. We present optimized algorithm in two platforms. Our GPU-based implementation exploits the parallelism of SIMD architectures to allow for maximum speed-up by a careful optimization of storage and data transfer, whereas our CPU-based implementation uses efficient sparse matrix computations and is suitable for large sparse data sets. For the community detection problem, we demonstrate accuracy and computational efficiency on Facebook, Yelp and DBLP data sets, and for the topic modeling problem, we also demonstrate good performance on the New York Times data set. We compare our results to the state-of-the-art algorithms such as the variational method, and report a gain of accuracy and a gain of several orders of magnitude in the execution time. We give efficient algorithms for volume sampling, i.e., for picking @math -subsets of the rows of any given matrix with probabilities proportional to the squared volumes of the simplices defined by them and the origin (or the squared volumes of the parallelepipeds defined by these subsets of rows). In other words, we can efficiently sample @math -subsets of @math with probabilities proportional to the corresponding @math by @math principal minors of any given @math by @math positive semi definite matrix. This solves an open problem from the monograph on spectral algorithms by Kannan and Vempala (see Section @math of KV , also implicit in BDM, DRVW ). Our first algorithm for volume sampling @math -subsets of rows from an @math -by- @math matrix runs in @math arithmetic operations (where @math is the exponent of matrix multiplication) and a second variant of it for @math -approximate volume sampling runs in @math arithmetic operations, which is almost linear in the size of the input (i.e., the number of entries) for small @math . Our efficient volume sampling algorithms imply the following results for low-rank matrix approximation: (1) Given @math , in @math arithmetic operations we can find @math of its rows such that projecting onto their span gives a @math -approximation to the matrix of rank @math closest to @math under the Frobenius norm. This improves the @math -approximation of Boutsidis, Drineas and Mahoney BDM and matches the lower bound shown in DRVW . The method of conditional expectations gives a algorithm with the same complexity. The running time can be improved to @math at the cost of losing an extra @math in the approximation factor. (2) The same rows and projection as in the previous point give a @math -approximation to the matrix of rank @math closest to @math under the spectral norm. In this paper, we show an almost matching lower bound of @math , even for @math . Strassen's matrix multiplication (MM) has benefits with respect to any (highly tuned) implementations of MM because Strassen's reduces the total number of operations. Strassen achieved this operation reduction by replacing computationally expensive MMs with matrix additions (MAs). For architectures with simple memory hierarchies, having fewer operations directly translates into an efficient utilization of the CPU and, thus, faster execution. However, for modern architectures with complex memory hierarchies, the operations introduced by the MAs have a limited in-cache data reuse and thus poor memory-hierarchy utilization, thereby overshadowing the (improved) CPU utilization, and making Strassen's algorithm (largely) useless on its own. In this paper, we investigate the interaction between Strassen's effective performance and the memory-hierarchy organization. We show how to exploit Strassen's full potential across different architectures. We present an easy-to-use adaptive algorithm that combines a novel implementation of Strassen's idea with the MM from automatically tuned linear algebra software (ATLAS) or GotoBLAS. An additional advantage of our algorithm is that it applies to any size and shape matrices and works equally well with row or column major layout. Our implementation consists of introducing a final step in the ATLAS GotoBLAS-installation process that estimates whether or not we can achieve any additional speedup using our Strassen's adaptation algorithm. Then we install our codes, validate our estimates, and determine the specific performance. We show that, by the right combination of Strassen's with ATLAS GotoBLAS, our approach achieves up to 30 22 speed-up versus ATLAS GotoBLAS alone on modern high-performance single processors. We consider and present the complexity and the numerical analysis of our algorithm, and, finally, we show performance for 17 (uniprocessor) systems. We study a statistical model for the tensor principal component analysis problem introduced by Montanari and Richard: Given a order- @math tensor @math of the form @math , where @math is a signal-to-noise ratio, @math is a unit vector, and @math is a random noise tensor, the goal is to recover the planted vector @math . For the case that @math has iid standard Gaussian entries, we give an efficient algorithm to recover @math whenever @math , and certify that the recovered vector is close to a maximum likelihood estimator, all with high probability over the random choice of @math . The previous best algorithms with provable guarantees required @math . In the regime @math , natural tensor-unfolding-based spectral relaxations for the underlying optimization problem break down (in the sense that their integrality gap is large). To go beyond this barrier, we use convex relaxations based on the sum-of-squares method. Our recovery algorithm proceeds by rounding a degree- @math sum-of-squares relaxations of the maximum-likelihood-estimation problem for the statistical model. To complement our algorithmic results, we show that degree- @math sum-of-squares relaxations break down for @math , which demonstrates that improving our current guarantees (by more than logarithmic factors) would require new techniques or might even be intractable. Finally, we show how to exploit additional problem structure in order to solve our sum-of-squares relaxations, up to some approximation, very efficiently. Our fastest algorithm runs in nearly-linear time using shifted (matrix) power iteration and has similar guarantees as above. The analysis of this algorithm also confirms a variant of a conjecture of Montanari and Richard about singular vectors of tensor unfoldings. This paper explores a seemingly counter-intuitive idea: the possibility of accelerating the solution of certain linear equations by adding even more equations to the problem. The basic insight is to trade-off problem size by problem structure. We test this idea on Toeplitz equations, in which case the expense of a larger set of equations easily leads to circulant structure. The idea leads to a very simple iterative algorithm, which works for a certain class of Toeplitz matrices, each iteration requiring only two circular convolutions. In the symmetric definite case, numerical experiments show that the method can compete with the preconditioned conjugate gradient method (PCG), which achieves O(nlogn) performance. Because the iteration does not converge for all Toeplitz matrices, we give necessary and sufficient conditions to ensure convergence (for not necessarily symmetric matrices), and suggest an efficient convergence test. In the positive definite case we determine the value of the free parameter of the circulant that leads to the fastest convergence, as well as the corresponding value for the spectral radius of the iteration matrix. Although the usefulness of the proposed iteration is limited in the case of ill-conditioned matrices, we believe that the results show that the problem size problem structure trade-off deserves further study. The massive parallelism of graphics processing units (GPUs) oers tremendous performance in many high-performance computing applications. While dense linear algebra readily maps to such platforms, harnessing this potential for sparse matrix computations presents additional challenges. Given its role in iterative methods for solving sparse linear systems and eigenvalue problems, sparse matrix-vector multiplication (SpMV) is of singular importance in sparse linear algebra. In this paper we discuss data structures and algorithms for SpMV that are eciently implemented on the CUDA platform for the ne-grained parallel architecture of the GPU. Given the memory-bound nature of SpMV, we emphasize memory bandwidth eciency and compact storage formats. We consider a broad spectrum of sparse matrices, from those that are well-structured and regular to highly irregular matrices with large imbalances in the distribution of nonzeros per matrix row. We develop methods to exploit several common forms of matrix structure while oering alternatives which accommodate greater irregularity. On structured, grid-based matrices we achieve performance of 36 GFLOP s in single precision and 16 GFLOP s in double precision on a GeForce GTX 280 GPU. For unstructured nite-element matrices, we observe performance in excess of 15 GFLOP s and 10 GFLOP s in single and double precision respectively. These results compare favorably to prior state-of-the-art studies of SpMV methods on conventional multicore processors. Our double precision SpMV performance is generally two and a half times that of a Cell BE with 8 SPEs and more than ten times greater than that of a quad-core Intel Clovertown system. We address the numerical problem of recovering large matrices of low rank when most of the entries are unknown. We exploit the geometry of the low-rank constraint to recast the problem as an unconstrained optimization problem on a single Grassmann manifold. We then apply second-order Riemannian trust-region methods (RTRMC 2) and Riemannian conjugate gradient methods (RCGMC) to solve it. A preconditioner for the Hessian is introduced that helps control the conditioning of the problem and we detail preconditioned versions of Riemannian optimization algorithms. The cost of each iteration is linear in the number of known entries. The proposed methods are competitive with state-of-the-art algorithms on a wide range of problem instances. In particular, they perform well on rectangular matrices. We further note that second-order and preconditioned methods are well suited to solve badly conditioned matrix completion tasks. Many problems in signal processing and statistical inference involve finding sparse solutions to under-determined, or ill-conditioned, linear systems of equations. A standard approach consists in minimizing an objective function which includes a quadratic (squared ) error term combined with a sparseness-inducing regularization term. Basis pursuit, the least absolute shrinkage and selection operator (LASSO), wavelet-based deconvolution, and compressed sensing are a few well-known examples of this approach. This paper proposes gradient projection (GP) algorithms for the bound-constrained quadratic programming (BCQP) formulation of these problems. We test variants of this approach that select the line search parameters in different ways, including techniques based on the Barzilai-Borwein method. Computational experiments show that these GP approaches perform well in a wide range of applications, often being significantly faster (in terms of computation time) than competing methods. Although the performance of GP methods tends to degrade as the regularization term is de-emphasized, we show how they can be embedded in a continuation scheme to recover their efficient practical performance. This article presents a fast solver for the dense "frontal" matrices that arise from the multifrontal sparse elimination process of 3D elliptic PDEs. The solver relies on the fact that these matrices can be efficiently represented as a hierarchically off-diagonal low-rank (HODLR) matrix. To construct the low-rank approximation of the off-diagonal blocks, we propose a new pseudo-skeleton scheme, the boundary distance low-rank approximation, that picks rows and columns based on the location of their corresponding vertices in the sparse matrix graph. We compare this new low-rank approximation method to the adaptive cross approximation (ACA) algorithm and show that it achieves better speedup specially for unstructured meshes. Using the HODLR direct solver as a preconditioner (with a low tolerance) to the GMRES iterative scheme, we can reach machine accuracy much faster than a conventional LU solver. Numerical benchmarks are provided for frontal matrices arising from 3D finite element problems corresponding to a wide range of applications. This paper describes a novel approach of packing sparse convolutional neural networks for their efficient systolic array implementations. By combining subsets of columns in the original filter matrix associated with a convolutional layer, we increase the utilization efficiency of the systolic array substantially (e.g., 4x) due to the increased density of nonzeros in the resulting packed filter matrix. In combining columns, for each row, all filter weights but one with the largest magnitude are pruned. We retrain the remaining weights to preserve high accuracy. We demonstrate that in mitigating data privacy concerns the retraining can be accomplished with only fractions of the original dataset (e.g., 10 for CIFAR-10). We study the effectiveness of this joint optimization for both high utilization and classification accuracy with ASIC and FPGA designs based on efficient bit-serial implementations of multiplier-accumulators. We present analysis and empirical evidence on the superior performance of our column combining approach against prior arts under metrics such as energy efficiency (3x) and inference latency (12x). Minimization with orthogonality constraints (e.g., (X^ X = I )) and or spherical constraints (e.g., ( x _2 = 1 )) has wide applications in polynomial optimization, combinatorial optimization, eigenvalue problems, sparse PCA, p-harmonic flows, 1-bit compressive sensing, matrix rank minimization, etc. These problems are difficult because the constraints are not only non-convex but numerically expensive to preserve during iterations. To deal with these difficulties, we apply the Cayley transform—a Crank-Nicolson-like update scheme—to preserve the constraints and based on it, develop curvilinear search algorithms with lower flops compared to those based on projections and geodesics. The efficiency of the proposed algorithms is demonstrated on a variety of test problems. In particular, for the maxcut problem, it exactly solves a decomposition formulation for the SDP relaxation. For polynomial optimization, nearest correlation matrix estimation and extreme eigenvalue problems, the proposed algorithms run very fast and return solutions no worse than those from their state-of-the-art algorithms. For the quadratic assignment problem, a gap 0.842 to the best known solution on the largest problem “tai256c” in QAPLIB can be reached in 5 min on a typical laptop. The problem of recovering a sparse signal x Rn from a relatively small number of its observations of the form y = Ax Rk, where A is a known matrix and k « n, has recently received a lot of attention under the rubric of compressed sensing (CS) and has applications in many areas of signal processing such as data cmpression, image processing, dimensionality reduction, etc. Recent work has established that if A is a random matrix with entries drawn independently from certain probability distributions then exact recovery of x from these observations can be guaranteed with high probability. In this paper, we show that Toeplitz-structured matrices with entries drawn independently from the same distributions are also sufficient to recover x from y with high probability, and we compare the performance of such matrices with that of fully independent and identically distributed ones. The use of Toeplitz matrices in CS applications has several potential advantages: (i) they require the generation of only O(n) independent random variables; (ii) multiplication with Toeplitz matrices can be efficiently implemented using fast Fourier transform, resulting in faster acquisition and reconstruction algorithms; and (iii) Toeplitz-structured matrices arise naturally in certain application areas such as system identification. In image processing, the Rudin-Osher-Fatemi (ROF) model [L. Rudin, S. Osher, and E. Fatemi, Phys. D, 60 (1992), pp. 259-268] based on total variation (TV) minimization has proven to be very useful. So far many researchers have contributed to designing fast numerical schemes and overcoming the nondifferentiability of the model. Methods considered to be particularly efficient for the ROF model include the Chan-Golub-Mulet (CGM) primal-dual method [T.F. Chan, G.H. Golub, and P. Mulet, SIAM J. Sci. Comput., 20 (1999), pp. 1964-1977], Chambolle's dual method [A. Chambolle, J. Math. Imaging Vis., 20 (2004), pp. 89-97], the splitting and quadratic penalty-based method [Y. Wang, J. Yang, W. Yin, and Y. Zhang, SIAM J. Imaging Sci., 1 (2008), pp. 248-272], and the split Bregman iteration [T. Goldstein and S. Osher, SIAM J. Imaging Sci., 2 (2009), pp. 323-343], as well as the augmented Lagrangian method [X.C. Tai and C. Wu, Lecture Notes in Comput. Sci. 5567, Springer-Verlag, Berlin, 2009, pp. 502-513]. In this paper, we first review the augmented Lagrangian method for the ROF model and then provide some convergence analysis and extensions to vectorial TV and high order models. All the algorithms and analysis will be presented in the discrete setting, which is much clearer for practical implementation than the continuous setting as in Tai and Wu, above. We also present, in the discrete setting, the connections between the augmented Lagrangian method, the dual methods, and the split Bregman iteration. Using our extensions and observations, we can easily figure out CGM and the split Bregman iteration for vectorial TV and high order models, which, to the best of our knowledge, have not been presented in the literature. Numerical examples demonstrate the efficiency and accuracy of our method, especially in the image deblurring case. In this paper we present a fast direct solver for certain classes of dense structured linear systems that works by first converting the given dense system to a larger system of block sparse equations and then uses standard sparse direct solvers. The kind of matrix structures that we consider are induced by numerical low rank in the off-diagonal blocks of the matrix and are related to the structures exploited by the fast multipole method (FMM) of Greengard and Rokhlin. The special structure that we exploit in this paper is captured by what we term the hierarchically semiseparable (HSS) representation of a matrix. Numerical experiments indicate that the method is probably backward stable. We describe new data structures for full and packed storage of dense symmetric triangular arrays that generalize both. Using the new data structures one is led to several new algorithms that save "half" the storage for symmetric matrices and outperform the current blocked based level 3 algorithms in LAPACK. We concentrate on the simplest forms of the new algorithms and show they are a direct generalization of LINPACK. This means that level 3 BLAS's are not required to obtain level 3 performance. The replacement for Level 3 BLAS are so-called kernel routines, see [1], and on IBM platforms they are producible from simple textbook type codes, by the XLF Fortran compiler. In the sequel I will label these "vanilla" codes. On Power3 with a peak performance of 800 MFlops, the results for Cholesky factorization at order n ? 200 is over 720 MFlops and then reaches 735 MFlops at n = 400. Using conventional full format LAPACK DPOTRF with ESSL BLAS's one first gets to 600 MFlops at n ? 600 and only reaches a peak of 620 MFlops. For this result we used simple square blocked full matrix data formats where the blocks themselves are stored in column major (Fortran) order or row major (C) order. The simple algorithms of LU factorization with partial pivoting for this new data format is a direct generalization of LINPACK algorithm DGEFA. Again, no conventional level 3 BLAS's are required; the replacements are again so-called kernel routines. Programming for squared blocked full matrix format can be accomplished in standard Fortran through the use of three and four dimensional arrays. Thus, no new compiler support is necessary. Also we mention that other more complicated algorithms are possible; e.g., recursive ones. The recursive algorithms are also easily programmed via the use of tables that address where the blocks are stored in the two dimensional recursive block array. Finally, we describe block hybrid formats. Doing so allows one to use no additional storage over conventional (full and packed) matrix storage. This means the new algorithms are completely portable. Many data analysis applications deal with large matrices and involve approximating the matrix using a small number of “components.” Typically, these components are linear combinations of the rows and columns of the matrix, and are thus difficult to interpret in terms of the original features of the input data. In this paper, we propose and study matrix approximations that are explicitly expressed in terms of a small number of columns and or rows of the data matrix, and thereby more amenable to interpretation in terms of the original data. Our main algorithmic results are two randomized algorithms which take as input an @math matrix @math and a rank parameter @math . In our first algorithm, @math is chosen, and we let @math , where @math is the Moore-Penrose generalized inverse of @math . In our second algorithm @math , @math , @math are chosen, and we let @math . ( @math and @math are matrices that consist of actual columns and rows, respectively, of @math , and @math is a generalized inverse of their intersection.) For each algorithm, we show that with probability at least @math , @math , where @math is the “best” rank- @math approximation provided by truncating the SVD of @math , and where @math is the Frobenius norm of the matrix @math . The number of columns of @math and rows of @math is a low-degree polynomial in @math , @math , and @math . Both the Numerical Linear Algebra community and the Theoretical Computer Science community have studied variants of these matrix decompositions over the last ten years. However, our two algorithms are the first polynomial time algorithms for such low-rank matrix approximations that come with relative-error guarantees; previously, in some cases, it was not even known whether such matrix decompositions exist. Both of our algorithms are simple and they take time of the order needed to approximately compute the top @math singular vectors of @math . The technical crux of our analysis is a novel, intuitive sampling method we introduce in this paper called “subspace sampling.” In subspace sampling, the sampling probabilities depend on the Euclidean norms of the rows of the top singular vectors. This allows us to obtain provable relative-error guarantees by deconvoluting “subspace” information and “size-of- @math ” information in the input matrix. This technique is likely to be useful for other matrix approximation and data analysis problems.