
Matrix decompositions are fundamental tools in the area of applied
mathematics, statistical computing, and machine learning. In particular,
lowrank 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 highdimensional data
matrix. The smaller matrix is then used to compute the desired lowrank
approximation. Such algorithms are shown to be computationally efficient for
approximating matrices with lowrank 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 multiscale
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 nonnegativity, 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
multicore 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 nonsmooth minimization problem directly. In contrast,
the approximations we propose allow us to replace the generalized nonsmooth
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 rankrevealing
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 GramSchmidt 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 columnpivoted QR, and is sometimes much
better. Computational speed is also improved substantially, in particular on
GPU architectures.