Randomized GPU Algorithms for the Construction of Hierarchical Matrices from Matrix-Vector Operations

Abstract

Randomized algorithms for the generation of low rank approximations of large dense matrices have become popular methods in scientific computing and machine learning. In this paper, we extend the scope of these methods and present batched GPU randomized algorithms for the efficient generation of low rank representations of large sets of small dense matrices, as well as their generalization to the construction of hierarchically low rank symmetric H 2 matrices with general partitioning structures. In both cases, the algorithms need to access the matrices only through matrix-vector multiplication operations which can be done in blocks to increase the arithmetic intensity and substantially boost the resulting performance. The batched GPU kernels are adaptive, allow nonuniform sizes in the matrices of the batch, and are more effective than SVD factorizations on matrices with fast decaying spectra. The hierarchical matrix generation consists of two phases, interleaved at every level of the matrix hierarchy. A first phase adaptively generates low rank approximations of matrix blocks through randomized matrix-vector sampling. A second phase accumulates and compresses these blocks into a hierarchical matrix that is incrementally constructed. The accumulation expresses the low rank blocks of a given level as a set of local low rank updates that are performed simultaneously on the whole matrix allowing high-performance batched kernels to be used in the compression operations. When the ranks of the blocks generated in the first phase are too large to be processed in a single operation, the low rank updates can be split into smaller-sized updates and applied in sequence. Assuming representative rank k, the resulting matrix has optimal O(kN) asymptotic storage complexity because of the nested bases it uses. The ability to generate an H 2 matrix from matrix-vector products allows us to support a general randomized matrix-matrix multiplication operation, an important kernel in hierarchical matrix computations. Numerical experiments demonstrate the high performance of the algorithms and their effectiveness in generating hierarchical matrices to a desired target accuracy.

Learning Elliptic Partial Differential Equations with Randomized Linear Algebra

“…We start by defining a GP on the domain Y . Let R Y ×Y K be the restriction 5 of the covariance kernel K to the domain Y × Y , which is a continuous symmetric positive definite kernel so that GP(0, R Y ×Y K ) defines a GP on Y . We choose a target rank k ≥ 1, an oversampling parameter p ≥ 2, and form a quasimatrix…”

Hierarchical algorithms on hierarchical architectures

“…Performance is obtained because the large amount of compute-intensive factorizations, both QR and SVD, that are performed at every level can be efficiently executed by batched kernels. We have developed batched QR and batched adaptive randomized SVD operations for this purpose [26,27].…”

“…We build on the hierarchical matrix vector and the low-rank update operations, to develop an algorithm for constructing a hierarchical matrix approximation of a 'black box' matrix that is accessible only via matrix-vector products [27]. The algorithm generalizes the popular randomized algorithms for generating low-rank approximations of large dense matrices [13] to the case of general hierarchical matrices H 2 .…”

Hierarchical matrix approximations for space-fractional diffusion equations

“…which produces the sum of a hierarchical matrix A of size N × N and a globally low rank matrix whose X and Y factors are of size N × k with k N . This operation can be efficiently implemented [30,54] by first adding the contributions of XY T to the various blocks of A at all levels, and recompressing the resulting sum algebraically as described earlier. The low rank update operation is a key routine for an operation that generates an explicit hierarchical matrix representation of an operator accessible only via matrix vector products.…”

“…This can be performed efficiently using Horner's method for evaluating polynomials. The resulting products are used in the level by level construction of the hierarchical matrix X p+1 [54]. We use methods of order 16 to construct the approximate inverse preconditioners used in the results shown below.…”