Ritesh Khan, Anderin Carson
7 min
This paper by Ritesh Khan and Anderin Carson introduces H_h-matrices, a clever evolution of hierarchical matrices (H-matrices) designed to tackle massive dense matrices from computational science—like those from PDEs, boundary element methods, or kernel methods in machine learning. Traditional dense matrices scale as O(N²) in storage and time, which is disastrous for large N. H-matrices approximate these using a tree structure (quad-trees in 2D, oct-trees in 3D) where 'admissible' off-diagonal blocks are replaced by low-rank approximations, slashing costs to near-linear while preserving accuracy.
Standard H-matrices use 'standard admissibility': a block is admissible (low-rank) if clusters are well-separated geometrically (dist(cluster centers) > diam(clusters)). This works great for smooth kernels but fails for near-field interactions, leaving dense blocks scattered. Weak admissibility (used in HODLR/HSS) is more aggressive, allowing low-rank everywhere off-diagonal—but it assumes unrealistic structure.
The innovation: hybrid admissibility. Use standard at coarse levels (large blocks, where separation is clear) and weak at fine levels (small blocks near diagonal). This confines dense blocks strictly to the diagonal, minimizing them. Proofs show H_h-matrices store less than standard H-matrices under mild conditions. Special cases recover pure H-matrices or HODLR seamlessly. Intuition: coarse levels capture global structure reliably; fine levels exploit local low-rank even without separation.
Modern hardware (NVIDIA Tensor Cores, TPUs) loves low-precision (fp16/bfloat16) for speed/memory wins. But naive low-precision kills accuracy in H-matrices due to rounding in low-rank factors.
The paper's rounding error analysis fills a gap: for admissible blocks, geometric clustering bounds the numerical rank, ensuring low-precision reps don't amplify errors. They derive explicit bounds for H_h-matrix-vector products (MVPs), proving stability if block precisions respect a threshold based on singular values and tree level.
Adaptive rule: For each block, compute min precision (e.g., fp16/fp32/fp64) via σ_min / ε_machine, where σ_min is the smallest kept singular value. Build the H_h-matrix by dynamically quantizing blocks—fp64 for tricky near-diagonal, fp16 for smooth far-off. Result: no accuracy loss, but huge savings.
Experiments on 2D/3D kernel matrices (e.g., Laplacian, Gaussian) show 11× storage reduction vs. uniform double-precision H-matrices. MVPs match reference accuracy (10^{-12} rel. error) at 1/10th memory. H_h structure enables this: fewer dense blocks mean more low-rank blocks ripe for quantization.
Why care? Exascale computing demands this. H-matrices already power fast solvers, preconditioners, N-body sims. Mixed-precision H_h-matrices accelerate them further on GPUs, enabling larger problems in climate modeling, quantum chemistry, or ML kernel approximations—without rework.
Takeaway: Hybrid admissibility + adaptive precision = practical path to memory-efficient, stable fast matrix ops. Bridges theory (error bounds) and hardware (low-precision accel).
Hierarchical matrices are data-sparse approximations of dense matrices that are widely used for fast matrix computations. Hierarchical matrices are built using a tree data structure, with low-rank blocks identified by various admissibility conditions, such as standard admissibility and weak admissibility. This paper introduces a novel hierarchical matrix framework, namely $\mathcal{H}_h$, based on a hybrid admissibility condition: we use the standard admissibility at the coarser levels (larger blocks) and the weak admissibility at the finer levels (smaller blocks). This hybrid strategy confines dense blocks only along the diagonal. We provide a criterion that ensures lower storage cost for $\mathcal{H}_h$-matrices compared to $\mathcal{H}$-matrices under the standard admissibility condition. We carry out a rounding error analysis of $\mathcal{H}_h$-matrices and show that the admissible blocks of $\mathcal{H}_h$-matrices can be represented in low precision (precision lower than the working precision) without degrading the overall approximation quality. We provide an explicit rule for dynamically selecting the precision of a given admissible block, thereby proposing an adaptive mixed precision algorithm for constructing and storing $\mathcal{H}_h$-matrices. Furthermore, we show that the use of mixed precision does not compromise the numerical stability and accuracy of the resulting $\mathcal{H}_h$-matrix-vector product. We perform a range of numerical experiments to validate our theoretical findings. Our numerical results show that the proposed adaptive mixed precision $\mathcal{H}_h$-matrices achieve significant storage reductions (up to $11 \times$) compared with uniform double precision standard admissibility-based $\mathcal{H}$-matrices, without compromising accuracy.
Alex: Right, so the strict one saves on compression costs for far groups but pays in extra full blocks nearby, while the loose one flips that?
Sam: Exactly. In a test with 262,000 points in 2D, the loose HODLR needed about twice the low-rank storage of the strict standard version for similar accuracy, because its blocks had higher numerical ranks—meaning more details to capture. The hybrid picks the best of both: strict rules high up where big blocks compress easily with few numbers, then loose rules low down where small blocks have low enough detail to compress well.
Alex: Huh, so the switch level is key to tipping storage in hybrid's favor... How do they decide where to switch for the best savings?
Sam: They propose a simple algorithm to find the optimal switch by building the standard version first, noting the extra neighbor storage at each level, then checking if switching earlier cuts more total space. It picks the level where the gain peaks, often one level before the end. In tests, this ensures the hybrid always uses less space than the standard setup.
Alex: And the experiments back that up across dimensions?
Sam: Yes—in 2D with over a quarter million points, the hybrid's total storage falls under the standard's for various accuracies. A 3D test with 64 thousand points confirms it, with clear gains up to about three and a half times in neighbor terms alone.
Alex: Okay, so hybrid nails the structure... But what's the logic behind blending that with adaptive precision per block to keep errors in check?
Sam: For each compressible block, they analyze rounding errors from lower-precision storage of its factors—like the skinny matrices that rebuild it. The rule sets bits per block proportional to its size, level, and relative strength in the whole matrix. This bounds the total error close to the target accuracy. Small or deep-tree blocks use fewer bits safely, since their mistakes add up less.
Alex: So no global accuracy loss, even mixed... How exactly does the system decide the bits for each one?
Sam: It first measures each block's relative size—its total "energy" compared to the whole matrix. Then, for a block at a certain level, it picks the lowest precision where the bit depth keeps errors under control. Their theorem ensures the sum of rounding mistakes across all blocks stays bounded by a small multiple of the target accuracy—no more than about four times in the bound.
Alex: And this holds for operations like multiplying the matrix by a vector?
Sam: Yes—they give an algorithm for that matrix-vector product, recursing down the tree and mixing precisions on the fly during computation. A figure with 6400 points shows it in action: for a tight target accuracy, most blocks drop to half-precision or below.
Alex: So the math checks out for storage and approximation... But how do they confirm it actually works in practice, like for real simulations?
Sam: They test on kernel matrices—grids where each entry depends only on the distance between pairs of points, like 1 over distance for simulating potentials around particles. Points are placed uniformly inside a cube or on a sphere surface, mimicking boundary element method setups in physics. Plots show adaptive mixed precision versions stay within predicted bounds for target accuracies, matching or beating uniform double-precision hybrids.
Alex: Okay, so errors controlled even when multiplying by a vector, like in iterative solvers. And storage tests shine?
Sam: Exactly—for 3D kernels at 125,000 points, adaptive hybrid uses about one-eleventh the space of standard double-precision versions at looser accuracies. The paper suggests this tackles memory crashes in large boundary simulations.
Alex: But to get those peak gains, they mostly set the switch level one step from the end—does hunting the true optimal one add much overhead?
Sam: It does—their simple algorithm builds a standard version first to scout levels, which takes extra passes through the tree. And it assumes a geometrically balanced tree, where splits keep clusters even-sized regardless of point distribution; uneven data might need adjustments.
Alex: Fair point—those assumptions keep the ranks low and predictable, but real-world point clouds aren't always uniform. Still, no more 100 gigabyte crashes for 125k-point 3D runs. A solid, practical step forward.
Sam: In sum, the hybrid matrices generalize prior approaches, their adaptive storage cuts costs substantially while preserving accuracy, and numerical results validate it across kernels and sizes. They plan to release the code publicly, which should aid further testing.
Alex: That's a meaningful contribution to handling massive simulations efficiently. Thanks, Sam—this has been a clear dive into making big math feasible on real hardware. Thanks for listening to ResearchPod.