PaperlayerView exact revision on arXiv

Source-linked AI summary

A distributed-memory package for dense Hierarchically Semi-Separable matrix computations using randomization

François-Henry Rouet, Xiaoye S. Li, Pieter Ghysels, Artem Napov

arXiv:1503.05464v2cs.MS

TL;DR

Dense matrices with low-rank off-diagonal blocks offer opportunities to reduce the cost of matrix operations, but parallel HSS software was scarce. This paper presents STRUMPACK, a distributed-memory package using randomized HSS compression and structured operations. Reported comparisons show STRUMPACK faster than Hsolver and ScaLAPACK, while some test matrices remain slower than ScaLAPACK.

  • Problem

    Dense matrices with low-rank off-diagonal blocks arise in applications including finite and boundary element methods, motivating lower-cost storage and matrix operations.

    Source: §1.1 · p. 1

  • Method

    STRUMPACK uses randomized sampling with adaptive compression to form HSS representations and provides distributed matrix-vector products, factorization, and solution routines.

    Source: §2.2 · p. 5, §5 · pp. 29–30

  • Results

    STRUMPACK was 2 to 6 times faster than Hsolver and 1.8 to 5.4 times faster than ScaLAPACK in the reported comparisons.

    Source: §4.6 · p. 27

  • Takeaways & Limitations

    The compression threshold lets STRUMPACK serve as a direct solver at full accuracy or as a robust preconditioner.

    Source: §5 · p. 29

    STRUMPACK was slower than ScaLAPACK on Schur100 and BEMAcoustics despite performing fewer floating-point operations.

    Source: §4.3 · p. 25

Abstract

from arXiv · show

We present a distributed-memory library for computations with dense structured matrices. A matrix is considered structured if its off-diagonal blocks can be approximated by a rank-deficient matrix with low numerical rank. Here, we use Hierarchically Semi-Separable representations (HSS). Such matrices appear in many applications, e.g., finite element methods, boundary element methods, etc. Exploiting this structure allows for fast solution of linear systems and/or fast computation of matrix-vector products, which are the two main building blocks of matrix computations. The compression algorithm that we use, that computes the HSS form of an input dense matrix, relies on randomized sampling with a novel adaptive sampling mechanism. We discuss the parallelization of this algorithm and also present the parallelization of structured matrix-vector product, structured factorization and solution routines. The efficiency of the approach is demonstrated on large problems from different academic and industrial applications, on up to 8,000 cores. This work is part of a more global effort, the STRUMPACK (STRUctured Matrices PACKage) software package for computations with sparse and dense structured matrices. Hence, although useful on their own right, the routines also represent a step in the direction of a distributed-memory sparse solver.

1 Introduction

The introduction motivates compressing low-rank off-diagonal blocks to reduce the costs of dense matrix operations and presents STRUMPACK as a distributed-memory HSS package for those computations.

1.1 Background

  • Low-rank off-diagonal blocks arise in finite and boundary element methods and other applications, and compressing them can reduce storage and matrix-operation costs.The passages connect these blocks in discretized PDE matrices to smooth Green’s functions.Source: §1.1 · p. 1
  • HSS and related low-rank representations were studied, but publicly available parallel libraries were scarce; prior tools included Hsolver and StruMF.Other cited packages include Hlib and MUMPS.Source: §1.1 · pp. 1–2

1.2 Contributions of this work

  • STRUMPACK provides distributed-memory HSS compression with randomized sampling, ULV-like linear-system solution, and HSS matrix-vector products.The introduction identifies these as core operations in the package.Source: §1.2 · p. 2
  • STRUMPACK supports any number of MPI processes and unbalanced hierarchical partitions, and is reported up to 6x faster than Hsolver.A flexible task-to-process mapping accommodates non-uniform partitions and unbalanced HSS trees.Source: §1.2 · p. 2
  • Adaptive parallel sampling supports black-box use when the matrix’s rank structure cannot be estimated beforehand.The authors describe this method as efficient and as increasing usability.Source: §1.2 · p. 2
  • The algorithms were evaluated on large problems from a wide range of academic and industrial applications using many cores.The paper also reports weak- and strong-scaling experiments.Source: §1.2 · p. 2

2 Background on Hierarchically Semi-Separable matrices

This section introduces HSS matrices, whose hierarchical low-rank off-diagonal structure supports compact representations and fast operations. It covers their nested-basis representation, randomized compression, matrix-vector products, and ULV-like factorization for solving linear systems, with costs depending on numerical rank.

2.1 Representation

  • HSS represents off-diagonal blocks with low-rank generators arranged on a cluster tree, with nested bases linking generators across levels.Leaf diagonal blocks remain uncompressed; low-rank off-diagonal blocks make the representation useful for reducing storage and enabling fast operations.Source: §2.1 · pp. 3–4

2.2 Compression with randomized sampling

  • Randomized compression samples the matrix using random vectors, then uses interpolative decomposition to construct HSS generators without requiring explicit access to every matrix entry.The method requires matrix-vector products and access to selected entries; it assumes a rank bound and uses d = r + p random-vector columns, with p a small oversampling parameter.Source: §2.2 · pp. 5–6

2.3 Matrix-vector product

  • After compression, HSS matrix-vector products cost O(rn), but compression costs O(rn^2), so products are amortized when performed repeatedly.The product algorithm traverses the HSS tree twice, accumulating V-generator actions and then applying U generators, coupling blocks, and diagonal blocks.Source: §2.3 · p. 7

2.4 ULV-like factorization

  • The ULV-like factorization solves HSS linear systems by exploiting generator structure derived from interpolative decomposition rather than using orthogonal transformations.The described approach eliminates O(n − r) unknowns through fast transformations and factors the remaining O(r) unknowns with standard LU in the original ULV formulation.Source: §2.4 · p. 7
  • Using more HSS tree levels avoids the one-stage O(n^3) cost: elimination leaves O(r) unknowns for a root-level LU factorization.The procedure applies transformations and LQ factorizations to diagonal blocks before the root-level LU.Source: §2.4 · p. 10
  • O(r^2n) is the stated complexity of the ULV-like factorization, which produces factors and intermediate matrices used in the solution phase.At non-root nodes, the algorithm outputs Qτ, Lτ, Wτ, and Ṽτ; the latter two support the solution phase.Source: §2.4 · p. 11

2.5 Solution using ULV factorization

  • O(rn) is the complexity of solving with the ULV factors, using a bottom-up traversal to transform the right-hand side and perform forward substitution.Updates from nodes are accumulated at the root; the intermediate solution is then transformed back to the original solution.Source: §2.5 · p. 12

3 Distributed-memory parallelism

This section develops and analyzes distributed-memory HSS algorithms, covering compression with adaptive sampling and its communication costs. It also describes tree- and node-level parallelism for factorization, solution, and matrix-vector products, including communication-cost comparisons for compression.

3.1 Task mapping

  • Static proportional mapping splits processes top-down across HSS-tree children, using interval size because node costs depend on unknown ranks.If ranks are unbalanced, compression workloads may be unbalanced; after compression, rank information can support remapping for later operations or related problems.Source: §3.1 · pp. 12–14

3.2 Parallel compression

  • Parallel compression forms samples with PBLAS PxGEMM, redistributes matrix blocks along a serialized postorder traversal, and processes serial subtrees and parallel nodes with suitable kernels.Serial subtrees use BLAS and LAPACK; parallel nodes use PBLAS and ScaLAPACK, including matrix products and Interpolative Decomposition.Source: §3.2 · pp. 14–15
  • In serial runs, modified xGEQP3 is typically two to three times faster than BLAS2 MGS, while the two options show little performance difference in parallel.The comparison is between modified xGEQP3 and the modified Gram-Schmidt implementation with column pivoting.Source: §3.2 · p. 16

3.3 Adaptive sampling mechanism

  • Adaptive sampling starts with few random vectors, adds more when a node's rank is too large, and resumes without recompressing completed work.In parallel, late notification lets branches progress while sufficient vectors remain; failed nodes are left untouched, and processes restart from the root after generating more vectors.Source: §3.3 · p. 16

3.4 Communication analysis

  • Communication analysis assumes uniform HSS rank, leaf diagonal blocks of size O(r), a complete binary tree, and a power-of-two process count.It counts messages and words on the critical path and analyzes sample computation, matrix distribution, and postorder generator construction.Source: §3.4 · p. 19
  • For latency, initial distribution dominates at small maximum rank, whereas tree traversal dominates at large rank.Source: §3.4 · p. 20
  • For bandwidth, tree traversal dominates matrix multiplication, which dominates initial distribution at large rank; the ordering reverses at small rank.Source: §3.4 · p. 20
  • For small rank, randomized compression communicates fewer messages and less volume than ScaLAPACK LU, while for large rank the comparison can reverse.Source: §3.4 · p. 20
  • Randomized compression has slightly higher latency complexity due to the log p distribution term and matrix-matrix-product latency.The authors propose reducing initial-distribution messages to O(p) and note potential benefits from communication-avoiding methods such as 2.5D matrix multiplication.Source: §3.4 · p. 20
  • Randomized compression communicates fewer words than non-randomized HSS compression because its sampling matrix-product and tree-traversal terms are smaller.Source: §3.4 · pp. 20–21

3.5 Parallel factorization, solution, and product

  • Factorization, triangular solution, and matrix-vector products use tree parallelism with proportional task mapping and node parallelism through PBLAS and ScaLAPACK.Serial subtrees are processed with sequential BLAS and LAPACK routines.Source: §3.5 · p. 21

4 Experimental results

Experiments across diverse structured matrices show how tree choice, compression tolerance, and adaptive sampling affect HSS computations. STRUMPACK achieves substantial speedups in several settings, while performance depends on matrix structure, rank, and communication costs.

4.1 Applications

  • STRUMPACK's experiments span Toeplitz, boundary-element, finite-difference, covariance, and H-matrix inputs, including cases with slowly growing or large maximum ranks.The test systems were run on NERSC's Hopper and Edison machines.Source: §4.1 · p. 21

4.2 General trees

  • A comb-shaped HSS tree reduces maximum rank to 700 and compresses much faster than a complete binary tree, whose unstructured block produces rank 10,000.Reweighting the process mapping to give right nodes 75% of their parent's processes reduces factorization time from 101.4 to 65.8 seconds, while compression time rises from 14.1 to 19.5 seconds.Source: §4.2 · pp. 22–23, Figure 6 · p. 22

4.3 Solving linear systems

  • For HMatrix, STRUMPACK is 72.9 times faster than ScaLAPACK at ε = 10^-4, but it is slower on Schur100 and BEMAcoustics.Speedups are 2.9x for BEMMultiSphere and 1.7x for Covar30; higher compression thresholds can speed compression and factorization but may reduce solution accuracy.Source: §4.3 · pp. 24–25, Table 3 · p. 24

4.4 Fast matrix-vector product

  • The QChemToeplitz experiment uses HSS matrix-vector products in a power method to compute the largest eigenvalue.The supplied passage describes the experiment's purpose but does not report its numerical outcome.Source: §4.4 · p. 25, Table 4 · p. 25

4.5 Adaptive-rank mechanism

  • Adaptive sampling's runtime effect depends on tree depth: it can cost time on a three-level tree yet run faster than fixed 3,000-vector sampling on a deeper tree.Starting with 1,000 vectors and adding 1,000 after each failure terminates after two restarts with 3,000 vectors; rank variation does not materially affect accuracy or representation size.Source: §4.5 · p. 26, Figure 8 · p. 26, Table 5 · p. 27

4.6 Scalability

  • STRUMPACK is 1.8–5.4 times faster than ScaLAPACK and 2–6 times faster than Hsolver across the Helmholtz frontal-matrix problems.The HSS factors are smaller than Hsolver's, and the maximum ranks are similar. In the strict weak-scaling experiment, compression time is reasonably constant when redistribution time is excluded; for the large-rank strong-scaling case, the gap over ScaLAPACK narrows as process counts rise.Source: §4.6 · pp. 27–28, Table 6 · p. 28, Figure 9 · p. 29

5 Conclusion

STRUMPACK compresses dense matrices into HSS form for operations such as solving linear systems and computing matrix-vector products. Its randomized compression uses adaptive sampling, and the package supports both direct-solver accuracy and robust-preconditioner use.

5 Conclusion

  • STRUMPACK uses randomized sampling with adaptive sampling to compress dense matrices into HSS form for solves and matrix-vector products.The adaptive mechanism supports black-box use when rank structure cannot be estimated a priori.Source: §5 · pp. 29–30
  • A compression threshold lets STRUMPACK serve either as a direct solver with full accuracy or as a robust preconditioner.Source: §5 · pp. 29–30
  • The package supports any number of MPI processes and different hierarchical matrix partitionings, and is open source under the BSD-LBNL license.The presented version is identified as STRUMPACK-Dense-0.9.0.Source: §5 · p. 30
  • Future work includes integrating the dense package into a sparse solver and adding matrix-free support based on matrix-vector products and selected-element access.The matrix-free feature is described as planned for a future STRUMPACK version.Source: §5 · p. 30

A Two-stage triangular solution process

The ULV solve applies transformations to form U^-1b, then performs forward triangular substitution bottom-up through the HSS hierarchy. Intermediate z_tau variables accumulate contributions from eliminated nodes before the final orthogonal transformation produces the solution.

A Two-stage triangular solution process

Loading 1503.05464v2…