REVIEW 1 major objections 1 minor 34 references
Adaptive Sketching Based Construction of H2 Matrices on GPUs
T0 review · 1 major / 1 minor · reviewed 2026-08-15 · deepseek-v4-flash
Pith's one-line read A bottom-up adaptive sketching algorithm constructs H2 matrices in linear time and memory, and runs on GPUs with up to 1000x speedups over existing methods.
desk verdict Strong GPU engineering for H2 construction, but the adaptive accuracy guarantee is heuristic and needs proof or explicit framing as such. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing object is the bottom-up skeletonization sweep together with the sample-propagation identity $Y^l = W^l Y^{l-1}$, where $W^l$ is a block unit-triangular matrix assembled from the interpolation matrices of interpolative decompositions at level $l$. This identity means that once the leaf level has subtracted its dense blocks, every higher level's sample matrix is obtained by simple block operations on the previous level's samples, so the black-box operator $K_{\mathrm{blk}}$ is invoked only for the initial sketch and for adaptive additions. The second essential piece is the adaptive convergence test: for each node, a QR factorization of the local residual sample matrix $Y_{\mathrm{loc}}^\tau$ is computed, and the node is declared converged when the smallest diagonal entry of the triangular factor falls below an absolute threshold $\epsilon_{\mathrm{abs}}$; samples are added in blocks until every node at the level passes. On the GPU side, the machinery is batching: all operations are marshaled into a few kernel launches per level, with the block-sparse row matrix multiplication split into at most $C_{\mathrm{sp}}$ batched non-uniform matrix multiplications to avoid atomics.
What would settle it
Run the adaptive algorithm on a matrix whose off-diagonal blocks need noticeably higher ranks at coarser levels than at the leaves — for example a three-dimensional integral-equation kernel with strong admissibility and a target tolerance of $10^{-6}$ — and record, at every node, whether the local QR-based test reports convergence while the final measured relative error $\|K_{\mathrm{comp}} - K\|/\|K\|$ exceeds the target; the claimed adaptive guarantee falls if the local test passes everywhere yet the global error is out of tolerance.
Extended reading notes
Core claim
The central claim, stated on the paper's own terms, is that the nested-basis $\mathcal{H}^2$ format can be built in linear complexity by a bottom-up sweep that never re-applies the black-box sketching operator after the initial samples are taken. At the leaf level the algorithm forms $Y_{\mathrm{loc}}^\tau = Y^\tau - \sum_{b\in \mathcal{N}_\tau} D_{\tau,b}\Omega_b$, removing the directly evaluated dense inadmissible blocks so that only admissible interactions remain in the sample matrix; a row interpolative decomposition of $Y_{\mathrm{loc}}^\tau$ then picks skeleton indices $\tilde{I}_\tau$ and defines the leaf basis. At every higher level, the samples and random vectors of a node are assembled from its children and the residual after subtracting the coupling-block contributions is again decomposed, yielding the transfer matrices $E_{\tau_1}, E_{\tau_2}$ and the skeleton set, with coupling matrices $B_{\tau,b} = K(\tilde{I}_\tau,\tilde{I}_b)$ evaluated directly. The key identity $Y^l = W^l Y^{l-1}$, with $W^l$ block unit-triangular, is what lets the sweep propagate samples without further sketching. The paper asserts that for a fixed tolerance the sample count remains $O(1)$ — 256 in the reported experiments — making the whole construction an $O(N)$ algorithm when the rank is $r = O(\log 1/\epsilon)$.
Load-bearing premise
The accuracy guarantee rests on the assumption that a node's approximation is converged once the smallest diagonal entry of the QR factor of its local sample matrix falls below an absolute threshold, and that checking this locally at every node, level by level, keeps the global relative error of the assembled $\mathcal{H}^2$ matrix within the target tolerance; the paper does not prove this and notes that its simple error compensation does not fully account for errors accumulated while sweeping up the tree.
Editorial extensions
If this is right
- For a fixed tolerance, construction costs $O(N)$ time and memory, so $\mathcal{H}^2$ representations of kernel and integral-equation matrices become practical at sizes where top-down sample counts would exhaust memory.
- The sample count stops growing with $N$ — 256 in the reported runs instead of thousands to tens of thousands — which removes the dominant cost of the black-box sketching operator and is the main source of the reported speedups.
- Because the construction is adaptive, users no longer need to guess the rank in advance; the algorithm adds sample blocks until the local convergence test is met.
- The same CPU/GPU code base with batched variable-size kernels can be used for low-rank updates of existing $\mathcal{H}^2$ matrices, a step that appears in multifrontal factorizations and Schur-complement updates.
Reading between the lines
- If the local QR convergence test is sound, the same adaptive criterion could be applied to other hierarchical formats (HODLR, HSS, fully black-box bottom-up skeletonization), potentially giving them the same sample-count reduction.
- A practitioner comparing methods may want to normalize by achieved accuracy: the baselines in the paper were constructed to $10^{-5}$ while the proposed algorithm targeted $10^{-6}$, so a time-to-accuracy comparison would put the speedup ratios in context.
- The batched variable-size marshaling pattern is reusable beyond $\mathcal{H}^2$: the same kernel design should accelerate other data-sparse formats and sparse multifrontal front updates on GPUs.
- The paper's admitted lack of a formal error analysis for the adaptive tolerance suggests a concrete research target: derive an a priori bound linking the per-node QR diagonal threshold to the global relative error, which would let users set $\epsilon_{\mathrm{abs}}$ with confidence.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper presents a bottom-up sketching-based algorithm for constructing H2 matrices with linear complexity, together with a GPU implementation based on batched dense linear algebra and batched entry evaluation. The algorithm requires a black-box matrix-vector product and an entry-evaluation routine, extends an earlier HSS construction method to strongly admissible H2 matrices, and includes an adaptive sampling mechanism that adds random vectors until a local QR-based convergence test is satisfied. Numerical experiments on covariance matrices, volume integral equation operators, and H2 update problems report speedups up to 13x over a CPU version, up to 1000x over a GPU top-down H2 construction in H2Opus, and up to 660x over ButterflyPACK, with memory scaling consistent with O(N). The paper claims to be the first GPU implementation of the bottom-up sketching-based H2 construction class.
Significance. If the correctness and complexity claims are fully substantiated, this is a valuable contribution: it demonstrates that bottom-up, strongly admissible H2 construction can be made adaptive, sample-efficient, and GPU-resident, and the batched implementation of non-uniform block-sparse operations is a nontrivial engineering achievement. The stated reduction in the number of random vectors relative to top-down methods is a concrete and useful advance. However, the central accuracy guarantee currently rests on an unproven local convergence test, and the headline speedup numbers are obtained from comparisons at different target tolerances; both issues need to be resolved before the performance claims can be accepted at face value.
major comments (1)
- [Section IV-B / Section V-C] The performance analysis in Section IV-B claims that the batched operations require only L=O(log N) kernel launches and that kernel-launch overhead is negligible, but the numerical breakdown in Fig. 7 shows that the convergence test (dominated by batched QR) takes up a significant fraction of GPU runtime, especially for smaller problems. The claim that the kernel-launch cost is negligible appears to be based on asymptotic counts rather than measured overhead, and the profiling data suggest that batched QR efficiency is a practical bottleneck. The authors should reconcile the asymptotic kernel-launch argument with the profiled time distribution, or restrict the claim accordingly.
minor comments (1)
- [Section I] The phrase 'It's also worth mentioning' should be changed to 'It is also worth mentioning' for formal style.
Circularity Check
No significant circularity: the H2 construction is an independent compression pipeline; deferred error analysis and a self-referential validation metric are caveats, not circular reductions.
full rationale
Algorithm 1 takes a black-box sketching operator K_blk(·) and an entry-evaluation function as inputs, and arrives at the H2 representation via interpolative decompositions of sketched sample blocks. No parameter is fitted to the reported accuracy, and no equation in the derivation is defined in terms of the output H2 matrix. The O(N) complexity claim is argued directly from the operation counts in Algorithm 1 (O(r^2 N) sampling and O(rN) entry evaluations under a fixed tolerance), not from a self-citation. The adaptive convergence test in Section III-B is a local QR-based criterion whose global error behavior is explicitly deferred to future work, and Section V-D admits that the error compensation 'does not fully account for the approximation errors as we sweep up the tree.' These are acknowledged limitations in the accuracy guarantee, but they are not circular: the algorithm does not define its reported error as the stopping criterion, and the experiments report measured errors below the target threshold. The relative error is measured against the provided sampler K_blk, and for covariance/IE experiments that sampler is itself an H2Opus H2 approximation. This makes the validation somewhat self-referential, but it is a stated experimental choice, not a load-bearing derivation; the construction is still an independent compression of the given oracle and is benchmarked against external packages (ButterflyPACK, STRUMPACK, H2Opus). Overall, no claimed result reduces by construction to its own inputs, so the paper is not circular.
Assumptions & free parameters
free parameters (3)
- Leaf size =
64, 128, 256 in experiments
- Initial sample block size =
256 (main runs), 32 (adaptive runs)
- Admissibility parameter eta =
0.7
assumptions (5)
- domain assumption The target matrix admits a strongly-admissible H2 representation with block ranks r = O(log 1/epsilon) and bounded sparsity constant C_sp.
- domain assumption The black-box operator K_blk computes Y = K Omega in O(Nd) time and entry evaluation is cheap.
- standard math Randomized sketching with Gaussian-like matrices provides a good column space with high probability, and column-pivoted QR gives a valid interpolative decomposition.
- ad hoc to paper The QR-diagonal convergence test at each node guarantees the global H2 approximation error stays below the target tolerance.
- domain assumption The matrix is symmetric and real-valued.
Cite this review
Pith. "Pith review of Adaptive Sketching Based Construction of H2 Matrices on GPUs." pith.science (2026). https://pith.science/paper/VXGKX7JZ
@misc{pith2026250616759,
author = {Pith},
title = {Pith review of: Adaptive Sketching Based Construction of H2 Matrices on GPUs},
year = {2026},
howpublished = {\url{https://pith.science/paper/VXGKX7JZ}},
note = {Machine review of arXiv:2506.16759}
}
abstract
We develop a novel linear-complexity bottom-up sketching-based algorithm for constructing a $H^2$ matrix, and present its high performance GPU implementation. The construction algorithm requires both a black-box sketching operator and an entry evaluation function. The novelty of our GPU approach centers around the design and implementation of the above two operations in batched mode on GPU with accommodation for variable-size data structures in a batch. The batch algorithms minimize the number of kernel launches and maximize the GPU throughput. When applied to covariance matrices, volume IE matrices and $H^2$ update operations, our proposed GPU implementation achieves up to $13\times$ speedup over our CPU implementation, and up to $1000\times$ speedup over an existing GPU implementation of the top-down sketching-based algorithm from the H2Opus library. It also achieves a $660\times$ speedup over an existing sketching-based $H$ construction algorithm from the ButterflyPACK library. Our work represents the first GPU implementation of the class of bottom-up sketching-based $H^2$ construction algorithms.
Figures
Figures from the paper (2 more)
Reference graph
Works this paper leans on
-
[1]
A sparse matrix arithmetic based onH-matrices. part i: Introduction to-matrices,
W. Hackbusch, “A sparse matrix arithmetic based onH-matrices. part i: Introduction to-matrices,”Computing, vol. 62, no. 2, pp. 89–108, 1999
work page 1999
-
[2]
M. Bebendorf and W. Hackbusch, “Existence ofH-matrix approximants to the inverse fe-matrix of elliptic operators withl ∞-coefficients,” Numerische Mathematik, vol. 95, pp. 1–28, 2003
work page 2003
-
[3]
Introduction to hierar- chical matrices with applications,
S. B ¨orm, L. Grasedyck, and W. Hackbusch, “Introduction to hierar- chical matrices with applications,”Engineering analysis with boundary elements, vol. 27, no. 5, pp. 405–422, 2003
2003
-
[4]
Y . Zhang, A. Gillman, and S. Veerapaneni, “A fast direct solver for integral equations on locally refined boundary discretizations and its application to multiphase flow simulations,”Advances in Computational Mathematics, vol. 48, no. 5, p. 63, 2022
work page 2022
-
[5]
Y . Tu, Z. Xu, and H. Yang, “Hierarchical interpolative factorization for self green’s function in 3d modified poisson-boltzmann equations,” Communications on Applied Mathematics and Computation, pp. 1–26, 2024
work page 2024
-
[6]
Randomized sparse direct solvers,
J. Xia, “Randomized sparse direct solvers,”SIAM Journal on Matrix Analysis and Applications, vol. 34, no. 1, pp. 197–227, 2013
work page 2013
-
[7]
P. Ghysels, S. L. Xiaoye, C. Gorman, and F.-H. Rouet, “A robust parallel preconditioner for indefinite systems using hierarchical matrices and randomized sampling,” in2017 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2017, pp. 897– 906
work page 2017
-
[8]
Scalable and memory-efficient kernel ridge regression,
G. Ch ´avez, Y . Liu, P. Ghysels, X. S. Li, and E. Rebrova, “Scalable and memory-efficient kernel ridge regression,” in2020 IEEE International parallel and distributed processing symposium (IPDPS). IEEE, 2020, pp. 956–965
work page 2020
Show all 34 references
-
[9]
Fast direct methods for Gaussian processes,
S. Ambikasaran, D. Foreman-Mackey, L. Greengard, D. W. Hogg, and M. O’Neil, “Fast direct methods for Gaussian processes,”IEEE transactions on pattern analysis and machine intelligence, vol. 38, no. 2, pp. 252–265, 2015
2015
-
[10]
A superfast algorithm for Toeplitz systems of linear equations,
S. Chandrasekaran, M. Gu, X. Sun, J. Xia, and J. Zhu, “A superfast algorithm for Toeplitz systems of linear equations,”SIAM Journal on Matrix Analysis and Applications, vol. 29, no. 4, pp. 1247–1266, 2008
2008
-
[11]
AnO(NlogN)fast direct solver for partial hierarchically semi-separable matrices: with application to radial basis function interpolation,
S. Ambikasaran and E. Darve, “AnO(NlogN)fast direct solver for partial hierarchically semi-separable matrices: with application to radial basis function interpolation,”Journal of Scientific Computing, vol. 57, pp. 477–501, 2013
2013
-
[12]
A fast solver for HSS representations via sparse matrices,
S. Chandrasekaran, P. Dewilde, M. Gu, W. Lyons, and T. Pals, “A fast solver for HSS representations via sparse matrices,”SIAM Journal on Matrix Analysis and Applications, vol. 29, no. 1, pp. 67–81, 2007
2007
-
[13]
A direct solver with O(N)complexity for integral equations on one-dimensional domains,
A. Gillman, P. M. Young, and P.-G. Martinsson, “A direct solver with O(N)complexity for integral equations on one-dimensional domains,” Frontiers of Mathematics in China, vol. 7, pp. 217–247, 2012
2012
-
[14]
The inverse fast multipole method: using a fast approximate direct solver as a preconditioner for dense linear systems,
P. Coulier, H. Pouransari, and E. Darve, “The inverse fast multipole method: using a fast approximate direct solver as a preconditioner for dense linear systems,”SIAM Journal on Scientific Computing, vol. 39, no. 3, pp. A761–A796, 2017
2017
-
[15]
Hierarchical interpolative factorization for elliptic operators: integral equations,
K. L. Ho and L. Ying, “Hierarchical interpolative factorization for elliptic operators: integral equations,”Communications on Pure and Applied Mathematics, vol. 69, no. 7, pp. 1314–1353, 2016
2016
-
[16]
H2pack: High-performance h2 matrix package for kernel matrices using the proxy point method,
H. Huang, X. Xing, and E. Chow, “H2pack: High-performance h2 matrix package for kernel matrices using the proxy point method,”ACM Trans. Math. Softw., vol. 47, no. 1, dec 2020
2020
-
[17]
H2Opus: a distributed-memory multi-GPU software package for non- local operators,
S. Zampini, W. Boukaram, G. Turkiyyah, O. Knio, and D. Keyes, “H2Opus: a distributed-memory multi-GPU software package for non- local operators,”Advances in Computational Mathematics, vol. 48, no. 3, p. 31, 2022
2022
-
[18]
DistributedH 2-matrices for boundary element methods,
S. B ¨orm, “DistributedH 2-matrices for boundary element methods,” ACM Trans. Math. Softw., vol. 49, no. 2, jun 2023
2023
-
[19]
ASKIT: An efficient, parallel library for high-dimensional kernel summations,
W. B. March, B. Xiao, C. D. Yu, and G. Biros, “ASKIT: An efficient, parallel library for high-dimensional kernel summations,”SIAM Journal on Scientific Computing, vol. 38, no. 5, pp. S720–S749, 2016
2016
-
[20]
Geometry-oblivious FMM for compressing dense SPD matrices,
C. D. Yu, J. Levitt, S. Reiz, and G. Biros, “Geometry-oblivious FMM for compressing dense SPD matrices,” inProceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2017, pp. 1–14
2017
-
[21]
Algorithmic patterns for H-matrices on many-core proces- sors,
P. Zaspel, “Algorithmic patterns for H-matrices on many-core proces- sors,”Journal of Scientific Computing, vol. 78, no. 2, pp. 1174–1206, 2019
2019
-
[22]
Fast construction of hierarchical matrix representation from matrix-vector multiplication,
L. Lin, J. Lu, and L. Ying, “Fast construction of hierarchical matrix representation from matrix-vector multiplication,”Journal of Computa- tional Physics, vol. 230, no. 10, pp. 4071–4087, 2011
2011
-
[23]
Randomized compression of rank- structured matrices accelerated with graph coloring,
J. Levitt and P.-G. Martinsson, “Randomized compression of rank- structured matrices accelerated with graph coloring,”arXiv preprint arXiv:2205.03406, 2022
2022 arXiv
-
[24]
Randomized strong recursive skele- tonization: Simultaneous compression and factorization ofH-matrices in the black-box setting,
A. Yesypenko and P. Martinsson, “Randomized strong recursive skele- tonization: Simultaneous compression and factorization ofH-matrices in the black-box setting,”arXiv preprint arXiv:2311.01451, 2023
2023
-
[25]
A recursive skeletoniza- tion factorization based on strong admissibility,
V . Minden, K. L. Ho, A. Damle, and L. Ying, “A recursive skeletoniza- tion factorization based on strong admissibility,”Multiscale Modeling & Simulation, vol. 15, no. 2, pp. 768–796, 2017
2017
-
[26]
Direct solution of generalH 2 -matrices with controlled accuracy and concurrent change of cluster bases for elec- tromagnetic analysis,
M. Ma and D. Jiao, “Direct solution of generalH 2 -matrices with controlled accuracy and concurrent change of cluster bases for elec- tromagnetic analysis,”IEEE Transactions on Microwave Theory and Techniques, vol. 67, no. 6, pp. 2114–2127, 2019
2019
-
[27]
AnO(N) distributed-memory parallel direct solver for planar integral equations,
T. Liang, C. Chen, P.-G. Martinsson, and G. Biros, “AnO(N) distributed-memory parallel direct solver for planar integral equations,” in2024 IEEE International Parallel and Distributed Processing Sympo- sium (IPDPS). IEEE, 2024, pp. 440–452
2024
-
[28]
Scalable linear time dense direct solver for 3-d problems without trailing sub-matrix dependencies,
Q. Ma, S. Deshmukh, and R. Yokota, “Scalable linear time dense direct solver for 3-d problems without trailing sub-matrix dependencies,” inSC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2022, pp. 1–12
2022
-
[29]
A fast randomized algorithm for computing a hier- archically semiseparable representation of a matrix,
P.-G. Martinsson, “A fast randomized algorithm for computing a hier- archically semiseparable representation of a matrix,”SIAM Journal on Matrix Analysis and Applications, vol. 32, no. 4, pp. 1251–1274, 2011
2011
-
[30]
Thrust: A productivity-oriented library for CUDA,
N. Bell and J. Hoberock, “Thrust: A productivity-oriented library for CUDA,” inGPU computing gems Jade edition. Elsevier, 2011, pp. 359–371
2011
-
[31]
KBLAS: An optimized library for dense matrix-vector multiplication on GPU accelerators,
A. Abdelfattah, D. Keyes, and H. Ltaief, “KBLAS: An optimized library for dense matrix-vector multiplication on GPU accelerators,” ACM Transactions on Mathematical Software (TOMS), vol. 42, no. 3, pp. 1–31, 2016
2016
-
[32]
Frame- work for batched and GPU-resident factorization algorithms to block Householder transformations,
A. Haidar, T. Dong, S. Tomov, P. Luszczek, and J. Dongarra, “Frame- work for batched and GPU-resident factorization algorithms to block Householder transformations,” inISC High Performance, Springer. Frankfurt, Germany: Springer, 07-2015 2015
2015
-
[33]
ButterflyPACK,
Y . Liu and USDOE, “ButterflyPACK,” 11 2018. [Online]. Available: https://www.osti.gov//servlets/purl/1564244
2018
-
[34]
STRUMPACK – STRUctured Matrices PACKage,
P. Ghysels and USDOE, “STRUMPACK – STRUctured Matrices PACKage,” 12 2014. [Online]. Available: https://www.osti.gov/biblio/ 1328126
2014
Reviewed August 15, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.