PaperlayerView exact revision on arXiv

Source-linked AI summary

Scaling H2 Factorization to Billion-Unknown Three-Dimensional Problems

Tianyu Liang, Chao Chen, James Demmel, Yang Liu

arXiv:2610.30131v1math.NAcs.NAmath-phmath.MP

TL;DR

Distributed strong recursive skeletonization faces substantial communication and parallelism challenges in three dimensions. The paper combines spatial coloring with delayed updates and streamed sketching in a distributed H2 factorization. The resulting configuration reduces factorization time and peak-process memory in the reported Laplace and Helmholtz comparisons.

  • Problem

    Distributed strong recursive skeletonization faces substantial communication costs and limited parallelism, especially for three-dimensional problems with larger stencils, box counts, and coarse-level ranks.

    Source: §I · pp. 1–2

  • Method

    The solver combines coordinate-parity spatial box coloring, on-demand Schur updates, streamed sketching, and BLAS-3 and shared-memory execution for distributed symmetric H2 factorization.

    Source: §I · p. 2, §IX · p. 9

  • Results

    6.9× for Laplace and 8.5× for Helmholtz factorization time reductions, with peak process memory reductions of 57% and 54%, respectively, versus the unsketched, eager configuration.

    Source: §IX · p. 10

  • Takeaways & Limitations

    The experiments demonstrate distributed strong recursive skeletonization for a 3D Laplace volume integral equation with 3.62 billion unknowns at billion-unknown scale.

    Source: §IX · p. 9

Abstract

from arXiv · show

Strong recursive skeletonization provides approximate direct solvers for dense matrices from integral equations and kernel methods, but distributed implementations face substantial communication costs and limited parallelism, especially in three dimensions. We present a distributed-memory $\mathcal{H}^2$ factorization solver for symmetric matrices that combines spatial box coloring, delayed Schur complement updates, streamed sketching-based interpolative decomposition, and BLAS-3 and shared-memory kernels. Spatial coloring enables every process to participate in each boundary stage, while delayed updates and streamed sketching reduce dense storage and communication. To our knowledge, this is the first distributed-memory strong recursive skeletonization factorization demonstrated for three-dimensional volume integral equations at billion-unknown scale. Experiments on up to 512 compute nodes demonstrate factorization of a three-dimensional Laplace system with 3.62 billion unknowns, alongside three-dimensional Helmholtz tests and a two-dimensional Mat\'ern system with 268 million unknowns. Combining streamed sketching and fully delayed updates reduces factorization time by $6.9\times$ for Laplace and $8.5\times$ for Helmholtz relative to the unsketched, eager configuration, while reducing maximum per-process resident memory by 57\% and 54\%, respectively.

I. INTRODUCTION

Dense-system factorization costs motivate rank-structured methods, but extending distributed strong recursive skeletonization to 3D raises storage, communication, and scheduling challenges. The paper introduces spatial coloring and storage-reducing optimizations, and evaluates them at billion-unknown scale.

I. INTRODUCTION

  • Dense factorization scales quadratically in storage for N unknowns, motivating rank-structured solvers that reduce costs to quasi-linear or linear complexity.Such factorizations can serve as approximate direct solvers for multiple right-hand sides or as preconditioners.Source: §I · p. 1

II. BACKGROUND

  • Strong recursive skeletonization alternates low-rank compression and elimination, passing retained skeleton degrees of freedom to coarser levels; bounded numerical ranks can yield linear complexity.The paper studies strongly admissible H2 factorizations of symmetric matrices, including real- and complex-symmetric kernels, without assuming positive definiteness.Source: §I · p. 1

I. INTRODUCTION

  • Distributed parallel efforts include shared-memory solvers, H2-ULV methods, iterative hierarchical solvers, and a 2D distributed RS-S solver that idles about three quarters of processes per boundary stage.The 2D solver partitions the domain among processes and separates interior from boundary elimination.Source: §I · p. 1
  • 3D stencils increase Schur-update targets and compression data, while growing coarse-level ranks and sparse boundary concurrency make storage, communication, and scheduling harder.A direct 3D process-grid coloring would leave seven of eight processes idle during each boundary stage; comparisons with 2D ranks are qualified by differing tolerances and experimental settings.Source: §I · pp. 1–2
  • The proposed solver colors spatial boxes so processes can contribute eligible boundary boxes at every stage, delays Schur updates, streams sketching, and uses hybrid parallelism.The evaluation includes a 3D Laplace system with 3,623,878,656 unknowns on 512 nodes, 3D Helmholtz tests, and a 2D Matérn system with 268,435,456 unknowns.Source: §I · p. 2

A. Hierarchical Domain Decomposition

The method organizes points in a spatial tree and uses interpolative compression to separate retained skeleton indices from redundant indices. Eliminating the redundant indices updates nearby interactions through Schur complements.

A. Hierarchical Domain Decomposition

  • A perfect spatial tree recursively subdivides a uniform point distribution into equal-sized children, with a bounded number of points per leaf; coordinate-parity coloring is restricted to this setting.The tree is a quadtree in 2D or octree in 3D, and the focus is 3D point sets.Source: A. Hierarchical Domain Decomposition · p. 2
  • Each box has at most 26 near-field neighbors and 98 distance-2 boxes in 3D, compared with 8 and 16 in 2D.Near-field boxes share a face, edge, or corner; the distance-2 far field comprises boxes at grid distance two.Source: A. Hierarchical Domain Decomposition · p. 2

B. Strong Skeletonization of a Box

  • For suitable smooth kernels and fixed separation geometry, off-diagonal interactions have bounded numerical rank at prescribed accuracy, enabling an interpolative decomposition to select skeleton and redundant indices.A rank-revealing QR computes the ID at O(mnk) cost for rank k.Source: B. Strong Skeletonization of a Box · p. 2
  • Compression applies an ID to current distance-2 interactions together with proxy-surface evaluations that capture the far-field column space.The current interactions may already include Schur-complement updates from previously eliminated boxes.Source: B. Strong Skeletonization of a Box · pp. 2–3
  • Eliminating redundant indices after ID decoupling leaves them interacting only with skeleton and near-field indices, then produces Schur-complement updates on those remaining variables.The redundant pivot block is assumed nonsingular and is factored with a pivoted factorization.Source: B. Strong Skeletonization of a Box · p. 3

C. Multi-Level Factorization

The factorization proceeds from leaves toward the root, repeatedly compressing and eliminating box variables before factoring the residual dense matrix. Applying the resulting inverse uses an upward transformation pass, pivot solves, and a reverse downward pass.

C. Multi-Level Factorization

  • At each level from leaves through level 1, child skeletons become the active variables for their parent boxes, and the residual matrix is assembled and factored at the root.The root residual uses a dense pivoted factorization.Source: C. Multi-Level Factorization · p. 3
  • The approximate factorization combines box transformations with a block-diagonal matrix of redundant pivot blocks and the root residual.The total transformation is the product of box transformations in their application order.Source: C. Multi-Level Factorization · p. 3
  • Applying the approximate inverse uses an upward pass, solves with redundant pivot and dense root factors, then a downward pass applying transposed transformations in reverse elimination order.The formulation assumes symmetry and exploits it in compression, factor storage, and update computation; nonsymmetric matrices require distinct left and right transformations.Source: C. Multi-Level Factorization · p. 3

D. Dependency Structure of H2 Factorization

The factorization's data dependencies require boxes to receive prior-neighbor contributions before compression or elimination, with updates confined to distance two. The distributed schedule uses spatial coloring for boundary elimination and exchanges data between color stages.

D. Dependency Structure of H2 Factorization

  • The sequential procedure builds the spatial tree, processes boxes from leaves upward, assembles the next-level active matrix, and factors the dense root residual.For each box it computes an ID, forms decoupled blocks, factors the redundant pivot block, and applies Schur updates.Source: D. Dependency Structure of H2… · p. 4
  • A box must receive contributions from previously eliminated neighbors before compression or elimination, and elimination changes interactions only within distance two.Compression reads current distance-2 blocks and proxy interactions; elimination writes Schur updates on skeleton and near-field indices.Source: D. Dependency Structure of H2… · pp. 3–4

III. DISTRIBUTED ALGORITHM

  • Coordinate-parity coloring schedules boundary elimination so every process with eligible boxes can participate at each stage, unlike 3D process-grid coloring that would idle seven of eight processes.At each level, processes own regular contiguous subdomains of boxes.Source: §III · p. 4
  • Boundary boxes are processed in 2^d color sets, with exchanges after each boundary color; same-color boxes are nonadjacent, while shared-target updates are handled separately.Interior boxes are processed after the final boundary exchange, and parity batches also support shared-memory execution.Source: §III · p. 4

1) Graph Coloring:

The algorithm colors spatial boxes by coordinate parity, then uses the colors to schedule distributed boundary elimination and shared-memory interior work. It assembles surviving skeletons at each level and can reduce the active process count when too few boxes remain per process.

1) Graph Coloring:

  • Coordinate-parity coloring gives four colors in 2D and eight in 3D, with same-color boxes nonadjacent and independently compressible or eliminable.Their update supports can still overlap on shared neighbor pairs, so concurrent contributions need deferred accumulation.Source: 1) Graph Coloring · p. 4

2) Elimination Order:

  • Interior boxes use the same color order as boundary boxes, exposing shared-memory parallelism after boundary communication is complete.Same-color boxes are processed concurrently, with shared-target updates handled separately.Source: 2) Elimination Order · pp. 4–5, 1) Graph Coloring · p. 4
  • After each level, surviving child skeletons are assembled into the next-level matrix; the algorithm redistributes it with fewer processes when each active process would own too few boxes.The dense root matrix is factored at level 0.Source: 2) Elimination Order · pp. 4–5
  • Boundary boxes are eliminated color by color, with updates exchanged among neighboring processes before the next dependent color or interior phase.Interior boxes are then processed in color batches after the final boundary exchange.Source: 2) Elimination Order · pp. 4–5

3) Communication Pattern:

The communication schedule exchanges required boundary updates before dependent colors and interior processing. To reduce the storage and communication burden of 3D Schur updates, the method can retain generators and form products on demand.

3) Communication Pattern:

  • The final boundary exchange delivers remote blocks or update representations needed for interior compression, after which interior elimination needs no further communication.Interior boxes have local near-field support, so their newly generated updates do not cross process boundaries.Source: 3) Communication Pattern · p. 5

IV. ON-DEMAND SCHUR UPDATES

  • An eliminated box updates up to 325 off-diagonal neighbor pairs in 3D, compared with 28 in 2D, motivating deferred Schur updates to reduce communication and peak memory.The baseline eager implementation computes and buffers dense updates after each color.Source: §IV · p. 5, 3) Communication Pattern · p. 5

A. A generator representation for Schur updates

  • The on-demand representation communicates the redundant diagonal block D_E and generator G_E, forming update products when needed for elimination, ID compression, or parent assembly.The retained source representation can also include coupling data needed when a target box has already been skeletonized.Source: A. A generator representation for… · p. 5

B. Two evaluation policies

  • Eager mode stores modified near- and far-field blocks, whereas on-demand mode retains generators, removes persistent distance-2 products, and accumulates near-field updates into persistent blocks.On-demand far-field contributions can be consumed directly in sketch space during a later ID.Source: B. Two evaluation policies · p. 5

V. SKETCHING-ACCELERATED INTERPOLATIVE DECOMPOSITION

The ID compression target grows costly in 3D because it can include contributions from as many as 98 distance-2 boxes. The section presents projection-based compression and a streamed variant that sketches update contributions from generators without materializing them.

V. SKETCHING-ACCELERATED INTERPOLATIVE DECOMPOSITION

  • Direct RRQR of the compression target can be expensive in 3D because it may include contributions from as many as 98 distance-2 boxes.The target combines active-matrix interactions and proxy rows, and its size drives both arithmetic and storage costs.Source: §V · p. 6

A. Materialize-then-sketch

  • Materialize-then-sketch applies RRQR to a smaller random projection of the assembled target, reducing QR work while retaining full-target storage.The experiments use s_B = b_B and a sparse CountSketch with ζ = 4.Source: A. Materialize-then-sketch · p. 6

B. Stream-and-sketch

  • Stream-and-sketch sketches original interaction blocks and proxy rows directly, then sketches generator-based Schur updates without materializing those update matrices.Contributions from multiple sources are accumulated into Y_B as each generator arrives.Source: B. Stream-and-sketch · p. 6

VI. SHARED-MEMORY PARALLELISM

Shared-memory kernels use color batches and within-box parallelism to exploit available threads across hierarchy levels. Deferred target-pair updates avoid concurrent writes and provide a deterministic accumulation order.

VI. SHARED-MEMORY PARALLELISM

  • Color batches expose parallelism across nonadjacent boxes, while shared-memory work can use all node cores at fine levels and split boxes at coarse levels.The implementation also specifies race-free, deterministic Schur-update accumulation.Source: §VI · p. 6, A. Color Batches and Deferred Updates · p. 6

A. Color Batches and Deferred Updates

  • Each affected neighbor pair is assigned to one thread, which combines batch-source contributions in a fixed Morton order to avoid write races and make GEMM accumulation deterministic.Different neighbor pairs are parallelized across threads.Source: A. Color Batches and Deferred Updates · pp. 6–7

B. BLAS-3 Kernels and Parallelism within a Box

  • Native-layout generator and sketch operations expose Level-3 BLAS work and avoid the temporary matrix and two full passes required by a transpose-based solve.At coarse levels, cooperating tasks can split right-side solves, large GEMMs, and kernel evaluation into fixed-size chunks.Source: B. BLAS-3 Kernels and Parallelism… · p. 7

VII. COMMUNICATION AND MEMORY ANALYSIS

The analysis models per-level communication and interaction storage for spatial coloring, comparing eager Schur updates with on-demand evaluation. Under its stated assumptions, on-demand evaluation reduces modeled communication payload and stored interaction memory.

VII. COMMUNICATION AND MEMORY ANALYSIS

  • The communication model assumes all Schur updates from a boundary box are needed remotely and uses fixed skeleton size k and active block size bℓ.It also sets |N| = 3^d − 1 and |M| = 5^d − 3^d, with redundant size rℓ = bℓ − k.Source: §VII · p. 7

A. Communication cost

  • For a coarse 3D box with |N| = 26, bℓ ≈ 8k, and rℓ ≈ 7k, the mode-0 to mode-1 payload ratio is approximately 15.This estimate assumes all neighbor updates are needed remotely and does not predict a 15× reduction in total communication time, which also includes latency.Source: A. Communication cost · p. 7

B. Stored interaction memory

  • The stored-interaction model excludes diagonal blocks, retained factorization data, temporary workspaces, communication buffers, and runtime overhead.It therefore does not model total retained memory or per-process peak memory.Source: B. Stored interaction memory · p. 7
  • With d = 3, |N| = 26, |M| = 98, bℓ = 8k, and rℓ = 7k, on-demand evaluation yields a 42% reduction in modeled interaction storage.The per-process peak-memory reduction can differ from this modeled interaction-storage reduction.Source: B. Stored interaction memory · p. 7

VIII. EXPERIMENTS

The experiments benchmark spatial coloring on Laplace, Helmholtz, and Matérn kernels using regular grids in two or three dimensions. They report factorization-time components and matvec and relative-residual errors for Laplace and Helmholtz.

VIII. EXPERIMENTS

  • The benchmark uses 3D Laplace and Helmholtz kernels and a 2D Matérn kernel, with regular grids discretizing a unit cube or unit square.All three discretizations satisfy A = A⊤; Helmholtz is complex symmetric and uses non-conjugating transposes in the factorization.Source: §VIII · pp. 7–8
  • The Matérn experiments use length scale ℓ = 0.1 and diagonal regularization σ²_n = 10^-3.Reported factorization time is decomposed into computation and accumulated boundary data-exchange communication time.Source: §VIII · p. 8
  • The optimized 3D timing runs use one MPI process per node, n_leaf = 63, and reduction threshold n_min = 8.The weak-scaling series use 64 or 128 OpenMP workers per process, while strong-scaling runs vary the worker count.Source: §VIII · p. 8
  • Table II reports matvec error e_mv and relative residual error e_res for spatial coloring on Laplace and Helmholtz kernels.For Helmholtz, dω denotes the wave divisor.Source: Table II · p. 8

A. Accuracy and Efficacy

The accuracy tests report low errors for Laplace and Helmholtz, while ill-conditioned Matérn results show a limitation of proxy-surface compression and a benefit from PCG. The combined sketching and delayed-update configuration reduces factorization time and resident memory, with the largest gains in the 3D tests.

A. Accuracy and Efficacy

  • Laplace matvec errors stayed below 10^-3, while Helmholtz residuals reached 6.55 × 10^-7 at the largest size with ϵ = 10^-3.The reported Laplace residuals range from 3.30 × 10^-3; Helmholtz tests use a wave divisor chosen to permit low numerical ranks and well-conditioned A.Source: A. Accuracy and Efficacy · p. 8
  • For ill-conditioned 2D Matérn, proxy-surface compression misses target matvec accuracy and direct-solver solution error is large; preconditioned CG reduces solution error to ≤1.2 × 10^-8 across three runs.The matvec error decreases as the compression tolerance is tightened.Source: A. Accuracy and Efficacy · p. 8

B. Ablation Study of Streaming and On-Demand Updates

  • 8.5× Helmholtz and 6.9× Laplace factorization-time reductions; peak-process peak memory fell 54% and 57%, respectively, versus the unsketched, eager configuration.Helmholtz changed from 848.27 s to 99.70 s and 188.15 GiB to 85.65 GiB; Laplace changed from 714.37 s to 103.72 s and 123.24 GiB to 52.91 GiB.Source: B. Ablation Study of Streaming and… · p. 9
  • For Matérn, the combined configuration gives 1.39× lower factorization time and about 10% lower RSS.The Matérn test uses N = 163842 = 268,435,456 and ϵ = 10^-11.Source: B. Ablation Study of Streaming and… · p. 9

C. Parallel Performance

The section reports billion-unknown 3D factorization tests and performance across weak- and strong-scaling experiments. Spatial coloring, on-demand Schur updates, and streamed sketching address distributed communication and storage costs.

1) Weak scaling:

  • In weak scaling, communication and idling drive much of the runtime growth, while coarse levels have few boxes and large active sets that limit concurrency.The compared large-scale configurations use different per-process workloads and thread counts, so they are not direct performance comparisons; average skeleton size rises from 44 at the leaf level to 97–98 at L − ℓ = 3.Source: 1) Weak scaling · p. 9

2) Strong scaling:

  • 8.3×: increasing OpenMP workers from 1 to 64 reduces total factorization time from 860 s to 104 s for the fixed-size Laplace strong-scaling run.Leaf-level time falls from 461 s to 35 s, while coarser levels benefit less; Table VI breaks down per-level factorization and communication times.Source: 2) Strong scaling · p. 9, Table VI · p. 10

IX. CONCLUSION

  • Spatial box coloring lets every process with eligible boxes participate at each boundary stage, while on-demand Schur updates, streamed sketching, BLAS-3 solves, and task splitting target storage, communication, and local execution.The algorithm is a distributed-memory strong recursive skeletonization method for symmetric H2 matrices.Source: §IX · p. 9
  • The 3D Laplace test factorized 3.62 billion unknowns on 512 nodes, the first demonstrated billion-unknown-scale distributed-memory strong recursive skeletonization factorization for 3D volume integral equations.The experiments also tested 3D Helmholtz and a 2D Matérn system with 268 million unknowns.Source: §IX · p. 9, 1) Weak scaling · p. 9, Figure 3 · p. 10
  • 6.9× for Laplace and 8.5× for Helmholtz factorization-time reductions, with peak per-process memory reductions of 57% and 54%, respectively, versus the unsketched, eager configuration.Table IV covers factorization time, boundary data-exchange time, and maximum per-process resident memory across Helmholtz, Laplace, and Matérn results.Source: §IX · p. 10, Table IV · p. 10
  • Future work includes communication-avoiding parallelization, GPU acceleration, and extensions to more general kernels and unstructured point distributions.Source: §IX · p. 10
Loading 2610.30131v1…