A randomized least squares solver for terabyte-sized dense overdetermined systems

Highlights

• A scalable randomized least-squares solver based on the Blendenpik algorithm in a distributed-memory setting is implemented and evaluated.

• The randomized solver employs four different random projection transforms that enable our algorithm to scale up to three times in terms of matrix sizes and provide a speedup of up to 7.5 times over a state-of-the-art scalable least-squares solver for tera-scale matrices.

• The computational cost of various stages of the randomized solver is determined by how overdetermined the input matrix is.

• The Blendenpik solver demonstrates significant strong and weak scaling and excellent numerical stability in terms of relative error for increasing matrix sizes. The backward error is comparable to the backward error achieved by the state-of-the-art solver.

Abstract

We present a fast randomized least-squares solver for distributed-memory platforms. Our solver is based on the Blendenpik algorithm, but employs multiple random projection schemes to construct a sketch of the input matrix. These random projection sketching schemes, and in particular the use of the randomized Discrete Cosine Transform, enable our algorithm to scale the distributed memory vanilla implementation of Blendenpik to terabyte-sized matrices and provide up to ×7.5 speedup over a state-of-the-art scalable least-squares solver based on the classic QR algorithm. Experimental evaluations on terabyte scale matrices demonstrate excellent speedups on up to 16,384 cores on a Blue Gene/Q supercomputer.

Introduction

The explosive growth of data over the past 20 years in various domains, ranging from physics and biological sciences to economics and social sciences, has led to a need to perform efficient and scalable analysis on such massive datasets. One of the most widely and routinely used primitives in statistical data analysis is least-squares regression: given a matrix $A\in {ℝ}^{m\times n}$ and a vector $\text{b}\in {ℝ}^m$, we seek to compute

$${\text{x}}^{*}=arg\underset{\text{x}\in {ℝ}^n}{min}{∥A\text{x}-\text{b}∥}_2.$$

Several algorithms have been proposed to solve large-scale least-squares problems in various distributed and parallel environments [1], returning solutions whose accuracy is close to machine precision. Recent approaches for large-scale least-squares problems include communication-avoiding factorizations [2], which scale well for a variety of shared memory and distributed memory platforms [3] and are based on the classic QR decomposition algorithm, an O(mn²) algorithm (assuming m ≥ n).

Recent years have witnessed an explosion of research on so-called Randomized Numerical Linear Algebra [4] (or RandNLA for short) algorithms, which leverage the power of randomization in order to perform standard matrix computations. One of the core problems that have been extensively researched in this emerging field is the least-squares regression problem of Eq. (1). Sarlos [5] and Drineas et al. [6] introduced the first randomized algorithms for this problem. These algorithms are based on the application of the sub-sampled Randomized Hadamard Transform to the columns of the input matrix in order to create a least-squares problem of smaller size that can be then solved exactly and whose solution provably approximates the solution of the original problem with very high probability. This was followed by the work of Rokhlin and Tygert [7], who used a subsampled Randomized Fourier Transform to form a preconditioner and then used a standard iterative solver to solve the preconditioned problem. At the same time, Avron et al. [8] introduced Blendenpik, an algorithm and a software package which was the first practical implementation of a RandNLA dense least-squares solver that consistently and comprehensively outperformed state-of-the-art implementations of the traditional QR-based O(mn²) algorithms. Since then, there has been extensive research on RandNLA algorithms for regression problems; see Yang et al. [9] for a recent survey.

So far, most research on randomized least-squares regression algorithms has focused on the single processor setting, with two important exceptions. Meng et al. [10] introduced LSRN, a distributed memory algorithm for least-squares problems based on random Gaussian projections. While the algorithm is still an O(mn²) algorithm, the benefits of randomization are apparent with respect to both constants in the asymptotic analysis, as well as its much improved efficiency on parallel environments. Yang et al. [9] consider RandNLA in a MapReduce-like framework called Spark. This framework is less appropriate for supercomputers, as it is fails to take advantage of their hardware architecture.

In this work, we explore the behavior of Blendenpik-type algorithms in a distributed memory setting. We show that variants of Blendenpik that use various batchwise transformations to compute preconditioners lead to implementations that are not only faster than state-of-the-art implementations of baseline least-squares solvers, but are also able to scale to much larger matrix sizes. In particular, we show that a Blendenpik based algorithm can solve least-squares regression problems with terabyte-sized (and larger) input matrices. Our implementation and experiments were run on AMOS,¹ the high-performance Blue Gene/Q supercomputer system at Rensselaer. AMOS has five racks, 5120 nodes (81,920 cores), and 81,920 GB of main memory. AMOS has a peak performance of one PetaFLOP (10¹⁵ floating point operations per second), and a 5-D torus network with 2 GB/s of bandwidth per link and 512 GB/s to 1 TB/s of bisection network bandwidth per rack, depending on the torus network configuration. Due to runtime constraints imposed by the scheduling system for each partition of AMOS, we limited our experiments to partitions containing up to 1024 nodes (16,384 cores). The Blue Gene/Q architecture supports a hybrid communication framework that uses the MPI (Message Passing Interface) [11] standard for distributed communication and multithreading using OpenMP [12].

Our main contributions in this paper are: (i) implementation of and experimentation with the Blendenpik algorithm on distributed-memory platforms; (ii) implementation of four randomized sketching transforms in the context of the Blendenpik algorithm; (iii) a batchwise transformation scheme that scales a distributed vanilla implementation of the Randomized Discrete Cosine Transform in the context of the Blendenpik algorithm by up to three times in terms of matrix sizes, and provides a speedup of up to 7.5 times over a state-of-the-art scalable least-squares solver for tera-scale matrices; (iv) a detailed evaluation of four randomized sketching transforms in the context of the Blendenpik algorithm and their parameters on the BG/Q, using up to 16,384 cores.

The full source code of our batchwise Blendenpik implementation is available for download at https://github.com/cjiyer/libskylark/tree/batchwiseblendenpik. The rest of this paper is organized as follows. Section 2 describes the Blendenpik algorithm and the various stages of the algorithm in detail. Section 3 highlights the distributed Blendenpik implementation in the Blue Gene/Q, as well as scalability issues in our implementation, and describes an approach to overcome them. Section 4 first describes experiments to tune our Blue Gene/Q environment for our evaluations and then discusses the outcome of our experimental evaluations.

Notation. Let A, B, … denote matrices and let x, y, z, … denote column vectors. Given a vector $\text{x}\in {ℝ}^m$, let ${∥\text{x}∥}_2^2={\sum }_{i=1}^m{\text{x}}_i^2$ be (the square of) its Euclidean norm; given a matrix $A\in {ℝ}^{m\times n}$, let ${∥A∥}_F^2={\sum }_{i,j=1}^{m,n}{A}_{\mathrm{ij}}^2$ be (the square of) its Frobenius norm. Let σ₁ ≥ σ₁ ≥ σ₂ ⋯ ≥ σ_r > 0 be the nonzero singular values of A, where $r=\mathrm{rank}\left(A\right)$ is the rank of the matrix A. Then, the condition number of A is equal to κ(A) = σ₁/σ_r.