[Paper Review] Hierarchical Matrix Operations on GPUs: Matrix-Vector Multiplication and Compression

Wajih Boukaram, George Turkiyyah+1 more|arXiv (Cornell University)|Feb 5, 2019

Parallel Computing and Optimization Techniques22 references4 citations

TL;DR

This paper presents high-performance GPU implementations of hierarchical matrix-vector multiplication (HMV) and matrix compression for $Θ^2$-matrices, leveraging flattened tree representations and batched linear algebra to achieve near-optimal $\mathcal{O}(n)$ complexity. It demonstrates over 550 GB/s bandwidth and 850 GFLOPS/s sustained performance on a P100 GPU, enabling efficient large-scale dense matrix operations in spatial statistics and scientific computing.

ABSTRACT

Hierarchical matrices are space and time efficient representations of dense matrices that exploit the low rank structure of matrix blocks at different levels of granularity. The hierarchically low rank block partitioning produces representations that can be stored and operated on in near-linear complexity instead of the usual polynomial complexity of dense matrices. In this paper, we present high performance implementations of matrix vector multiplication and compression operations for the $\mathcal{H}^2$ variant of hierarchical matrices on GPUs. This variant exploits, in addition to the hierarchical block partitioning, hierarchical bases for the block representations and results in a scheme that requires only $O(n)$ storage and $O(n)$ complexity for the mat-vec and compression kernels. These two operations are at the core of algebraic operations for hierarchical matrices, the mat-vec being a ubiquitous operation in numerical algorithms while compression/recompression represents a key building block for other algebraic operations, which require periodic recompression during execution. The difficulties in developing efficient GPU algorithms come primarily from the irregular tree data structures that underlie the hierarchical representations, and the key to performance is to recast the computations on flattened trees in ways that allow batched linear algebra operations to be performed. This requires marshaling the irregularly laid out data in a way that allows them to be used by the batched routines. Marshaling operations only involve pointer arithmetic with no data movement and as a result have minimal overhead. Our numerical results on covariance matrices from 2D and 3D problems from spatial statistics show the high efficiency our routines achieve---over 550GB/s for the bandwidth-limited mat-vec and over 850GFLOPS/s in sustained performance for the compression on the P100 Pascal GPU.

Motivation & Objective

• Address the lack of high-performance GPU algorithms for core $Θ^2$-matrix operations, particularly matrix-vector multiplication and compression.
• Overcome the performance challenges posed by irregular tree data structures in hierarchical matrices on GPU-architectures.
• Enable efficient, memory-optimal computation of large dense matrices—common in spatial statistics and PDEs—on GPUs with limited global memory.
• Develop a flattened tree representation that allows efficient batching of linear algebra operations on irregular hierarchical data.
• Achieve near-linear complexity ($\mathcal{O}(n)$) for both HMV and compression, matching the theoretical efficiency of $Θ^2$-matrices on GPU hardware.

Proposed method

• Flatten the hierarchical tree structure of $Θ^2$-matrices into contiguous arrays to enable batched execution of linear algebra kernels.
• Use pointer arithmetic to marshal irregularly laid-out data into contiguous batches for GPU kernels, minimizing overhead.
• Implement matrix-vector multiplication using batched GEMV operations on flattened low-rank blocks and hierarchical basis structures.
• Perform matrix compression via batched QR and SVD operations on coupling matrices and basis nodes, with truncation to maintain accuracy.
• Leverage CUBLAS batched GEMM routines for orthogonalization and projection phases, optimizing for compute and memory bandwidth.
• Separate tree traversal logic from linear algebra computation to cleanly decouple data layout from algorithmic kernels.

Experimental results

Research questions

• RQ1 Can hierarchical matrix operations be efficiently mapped to GPU architectures despite their inherently recursive and irregular tree-based data structures?
• RQ2 How can the $\mathcal{O}(n)$ storage and computational complexity of $\mathcal{H}^2$-matrices be fully exploited on GPUs through algorithmic restructuring?
• RQ3 What performance gains can be achieved by flattening hierarchical tree representations and using batched linear algebra on GPUs for HMV and compression?
• RQ4 To what extent can the memory-bound HMV and compute-bound compression operations achieve high bandwidth and sustained FLOP rates on modern GPUs?
• RQ5 How effective is algebraic compression in reducing memory footprint and improving subsequent HMV performance for large-scale covariance matrices?

Key findings

• The matrix-vector multiplication (HMV) operation achieves over 550 GB/s bandwidth on the P100 GPU, surpassing the STREAM bandwidth benchmark for large 3D spatial statistics problems.
• The matrix compression operation sustains over 850 GFLOPS/s on the P100 GPU, with the orthogonalization phase achieving more than 2,000 GFLOPS/s.
• The basis generation and truncation phases of compression operate at over 600 GFLOPS/s, demonstrating high compute efficiency on the GPU.
• For a 1M×1M covariance matrix, HMV completes in under 29 ms, achieving 78% of the theoretical peak memory bandwidth on the P100.
• Compression of a 1M×1M matrix from an initial analytical representation to an optimal algebraic form completes in under 1.7 seconds on the P100 GPU.
• Memory savings from compression are significant—especially in the low-rank portions—and lead to measurable improvements in subsequent HMV runtime.