• Matrix decompositions are fundamental tools in the area of applied mathematics, statistical computing, and machine learning. In particular, low-rank matrix decompositions are vital, and widely used for data analysis, dimensionality reduction, and data compression. Massive datasets, however, pose a computational challenge for traditional algorithms, placing significant constraints on both memory and processing power. Recently, the powerful concept of randomness has been introduced as a strategy to ease the computational load. The essential idea of probabilistic algorithms is to employ some amount of randomness in order to derive a smaller matrix from a high-dimensional data matrix. The smaller matrix is then used to compute the desired low-rank approximation. Such algorithms are shown to be computationally efficient for approximating matrices with low-rank structure. We present the R package rsvd, and provide a tutorial introduction to randomized matrix decompositions. Specifically, randomized routines for the singular value decomposition, (robust) principal component analysis, interpolative decomposition, and CUR decomposition are discussed. Several examples demonstrate the routines, and show the computational advantage over other methods implemented in R.
  • The conjugate gradient method is a widely used algorithm for the numerical solution of a system of linear equations. It is particularly attractive because it allows one to take advantage of sparse matrices and produces (in case of infinite precision arithmetic) the exact solution after a finite number of iterations. It is thus well suited for many types of inverse problems. On the other hand, the method requires the computation of the gradient. Here difficulty can arise, since the functional of interest to the given inverse problem may not be differentiable. In this paper, we review two approaches to deal with this situation: iteratively reweighted least squares and convolution smoothing. We apply the methods to a more generalized, two parameter penalty functional. We show advantages of the proposed algorithms using examples from a geotomographical application and for synthetically constructed multi-scale reconstruction and regularization parameter estimation.
  • We present a new algorithm and the corresponding convergence analysis for the regularization of linear inverse problems with sparsity constraints, applied to a new generalized sparsity promoting functional. The algorithm is based on the idea of iteratively reweighted least squares, reducing the minimization at every iteration step to that of a functional including only $\ell_2$-norms. This amounts to smoothing of the absolute value function that appears in the generalized sparsity promoting penalty we consider, with the smoothing becoming iteratively less pronounced. We demonstrate that the sequence of iterates of our algorithm converges to a limit that minimizes the original functional.
  • The manuscript describes efficient algorithms for the computation of the CUR and ID decompositions. The methods used are based on simple modifications to the classical truncated pivoted QR decomposition, which means that highly optimized library codes can be utilized for implementation. For certain applications, further acceleration can be attained by incorporating techniques based on randomized projections. Numerical experiments demonstrate advantageous performance compared to existing techniques for computing CUR factorizations.
  • RSVDPACK is a library of functions for computing low rank approximations of matrices. The library includes functions for computing standard (partial) factorizations such as the Singular Value Decomposition (SVD), and also so called "structure preserving" factorizations such as the Interpolative Decomposition (ID) and the CUR decomposition. The ID and CUR factorizations pick subsets of the rows/columns of a matrix to use as bases for its row/column space. Such factorizations preserve properties of the matrix such as sparsity or non-negativity, are helpful in data interpretation, and require in certain contexts less memory than a partial SVD. The package implements highly efficient computational algorithms based on randomized sampling, as described and analyzed in [N. Halko, P.G. Martinsson, J. Tropp, "Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions," SIAM Review, 53(2), 2011], and subsequent papers. This manuscript presents some modifications to the basic algorithms that improve performance and ease of use. The library is written in C and supports both multi-core CPU and GPU architectures.
  • We introduce and compare new compression approaches to obtain regularized solutions of large linear systems which are commonly encountered in large scale inverse problems. We first describe how to approximate matrix vector operations with a large matrix through a sparser matrix with fewer nonzero elements, by borrowing from ideas used in wavelet image compression. Next, we describe and compare approaches based on the use of the low rank SVD, which can result in further size reductions. We describe how to obtain the approximate low rank SVD of the original matrix using the sparser wavelet compressed matrix. Some analytical results concerning the various methods are presented and the results of the proposed techniques are illustrated using both synthetic data and a very large linear system from a seismic tomography application, where we obtain significant compression gains with our methods, while still resolving the main features of the solutions.
  • We present new convolution based smooth approximations to the absolute value function and apply them to construct gradient based algorithms such as the nonlinear conjugate gradient scheme to obtain sparse, regularized solutions of linear systems $Ax = b$, a problem often tackled via iterative algorithms which attack the corresponding non-smooth minimization problem directly. In contrast, the approximations we propose allow us to replace the generalized non-smooth sparsity inducing functional by a smooth approximation of which we can readily compute gradients and Hessians. The resulting gradient based algorithms often yield a good estimate for the sought solution in few iterations and can either be used directly or to quickly warm start existing algorithms.
  • This manuscript describes a technique for computing partial rank-revealing factorizations, such as, e.g, a partial QR factorization or a partial singular value decomposition. The method takes as input a tolerance $\varepsilon$ and an $m\times n$ matrix $A$, and returns an approximate low rank factorization of $A$ that is accurate to within precision $\varepsilon$ in the Frobenius norm (or some other easily computed norm). The rank $k$ of the computed factorization (which is an output of the algorithm) is in all examples we examined very close to the theoretically optimal $\varepsilon$-rank. The proposed method is inspired by the Gram-Schmidt algorithm, and has the same $O(mnk)$ asymptotic flop count. However, the method relies on randomized sampling to avoid column pivoting, which allows it to be blocked, and hence accelerates practical computations by reducing communication. Numerical experiments demonstrate that the accuracy of the scheme is for every matrix that was tried at least as good as column-pivoted QR, and is sometimes much better. Computational speed is also improved substantially, in particular on GPU architectures.