{"id":"12c1f8b4-ef92-4c46-9d22-16e90eaf7b04","arxiv_id":"2506.16759","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"A bottom-up sketching algorithm constructs H2 matrices in linear time and is implemented on GPUs with batched kernels, reporting up to 1000x speedup over the H2Opus top-down GPU construction.","lead":"The paper presents a fast way to build H2 matrices, a compressed format for dense matrices that appear in physics and statistics, by repeatedly sampling the matrix instead of computing every entry. The authors also implement the method on GPUs and report large speedups over existing libraries, making the compressed format practical at larger scales.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Adaptive stopping criterion (Alg. 1 lines 11/29) certifies only local rank deficiency of Y_loc_tau, not a global H2 error bound, and Sec. V-D admits sweep-up errors are not fully accounted for.","rationale":"The reader's weakest assumption identifies the same load-bearing concern: the adaptive convergence test lacks a proof that local QR-diagonal thresholds control the global H2 approximation error. This is the least secure pillar of the central claim. The O(N) complexity claim is defensible under the stated fixed-rank assumption, and the GPU speedups are empirical and plausible, so they are not the main risk. The accuracy guarantee, by contrast, is both deferred and explicitly hedged in Section V-D, and it is essential to the advertised 'adaptive' construction at a target tolerance. The numerical experiments are consistent with the claim, so the appropriate response is to keep the verdict CONDITIONAL rather than to reject. The surrogate error measurement against the H2Opus 1e-5 oracle is a related but secondary issue; it weakens the empirical evidence for the 1e-6 accuracy claim, but the primary gap is the missing error analysis of the adaptive stopping rule. No code release is a reproducibility concern, not a correctness concern, and does not change the verdict.","tokens_in":15308,"tokens_out":8630,"duration_ms":98416,"concrete_test":"Analytical check: write the two-level error recursion. For a parent tau with children nu1 and nu2, compare Y_loc_tau as formed in Line 27 with the exact residual K_admissible Omega: the difference is a sum of child-level ID errors propagated through the computed B's. Derive whether the operator norm of this difference is bounded by C * eps_l, with C independent of L, N, and the number of adaptive rounds. If the bound contains an uncontrolled term such as ||E_child|| * eps_child or a C that grows with level, the QR-diagonal test cannot certify the global tolerance. Numerical check: for N=2^14 with the covariance kernel, use a dense reference matrix rather than the H2Opus 1e-5 oracle, run Algorithm 1 adaptively at eps=1e-6 for leaf sizes 64/128/256, and compute true ||K - K_H||/||K||; any result above 1e-6 confirms the concern.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section III-B defines convergence for a node tau by computing the QR of the local sample matrix Y_loc_tau and checking whether the smallest diagonal of R is below eps_abs (Algorithm 1, Lines 11 and 29). This test certifies that the sampled block-row is numerically rank-deficient, not that the assembled H2 approximation K_H satisfies ||K - K_H||/||K|| <= eps. The construction then treats the ID factors and transfer matrices as exact when forming the next level's Y_loc_tau: Eq. (7) and the block elimination in (5)-(6) require the B and T factors to reproduce K_l exactly, but at higher levels those factors come from previous approximate IDs. Errors from child-level compression therefore enter the parent residual, and the local QR diagonal does not bound their contribution. The paper explicitly leaves the error analysis to future work (end of Sec. III-A) and admits in Sec. V-D that 'the simple error compensation scheme [does] not fully account for the approximation errors as we sweep up the tree.' Since the adaptive sampler is the mechanism that nominally guarantees the advertised 1e-6 accuracy, this is a load-bearing gap: the algorithm can stop adding samples at a node while the global error is still above tolerance, silently. Reported accuracies in Table II are below 1e-6, so the concern is not a demonstrated failure, but the central accuracy claim is currently unproven.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":15603,"tokens_out":3332,"duration_ms":38042,"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":[{"comment":"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.","section":"Section IV-B / Section V-C"}],"minor_comments":[{"comment":"The phrase 'It's also worth mentioning' should be changed to 'It is also worth mentioning' for formal style.","section":"Section I"}],"recommendation":"major_revision","confidential_remarks":"The main obstacle to acceptance is the missing global error guarantee for the adaptive stopping rule and the fact that the reported accuracy is measured relative to an approximate sampler. The engineering contribution and the experimental scaling are genuinely useful, but the correctness claim and the fairness of the speedup comparisons need to be addressed before publication. I would be willing to review a revised version that provides at least a sketch of the global error analysis or a modified stopping criterion with a verifiable error certificate, and that repeats the reference comparisons at equal tolerances."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Colleague,\n\nThe real news in this paper is the GPU engineering. The authors extend Martinsson's HSS sketching to strongly-admissible H2 matrices with adaptive sampling, and implement it with batched variable-size kernels that handle the non-uniform block structure of H2 on a GPU. The speedups are large: up to 13x over their CPU version, three orders of magnitude over H2Opus's top-down GPU construction, and 660x over ButterflyPACK. The O(N) scaling in time and memory is backed by experiments up to N=524288. This is a genuine contribution and likely the first GPU implementation of bottom-up sketching-based H2 construction.\n\nThe soft spots are real but not fatal. The adaptive convergence test in Section III-B checks the smallest diagonal of the QR factor of the local sample matrix Y_loc_tau. That certifies local rank deficiency of the sample block-row, not a global bound on the H2 approximation error. The authors explicitly defer the error analysis to future work, and in Section V-D admit that their error compensation does not fully account for errors as they sweep up the tree. Since the adaptive sampler is what nominally guarantees the 1e-6 accuracy, this is a gap. It is not a demonstrated failure - Table II reports errors below 1e-6 - but the guarantee is unproven. A rigorous analysis of the adaptive stopping criterion, or at least a clear statement that it is a heuristic, is needed.\n\nTwo smaller issues. First, the reported error is measured against the provided sampler, not the exact matrix. In the low-rank update experiment, that means the error is relative to the old H2 matrix, not the exact updated matrix. Second, the baselines H2Opus and ButterflyPACK were constructed to 1e-5 while the proposed method targets 1e-6, which makes the speedup comparison not perfectly fair, though I would be surprised if the three orders of magnitude vanish under a matched tolerance. No code is released, which limits reproducibility, but the paper gives enough detail to reimplement.\n\nOverall, this is a serious and useful paper. The GPU work is substantial, the numerical evidence is consistent with the complexity claims, and the limitations are honestly acknowledged. I would send it to peer review and ask for the error analysis to be either supplied or explicitly scoped out. The paper will be of interest to anyone working on H2 construction, randomized sketching for rank-structured matrices, or GPU solvers for dense structured systems.","headline":"Strong GPU engineering for H2 construction, but the adaptive accuracy guarantee is heuristic and needs proof or explicit framing as such.","tokens_in":16113,"tokens_out":3076,"would_cite":true,"duration_ms":29028,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"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.","keywords":["H2 matrix","adaptive sketching","bottom-up construction","GPU batched kernels","interpolative decomposition","nested basis","linear complexity","randomized numerical linear algebra"],"falsifier":"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.","tokens_in":15114,"feed_emoji":"⚡","tokens_out":13718,"duration_ms":121365,"temperature":0.7,"pith_summary":"The paper sets out to show that an $\\mathcal{H}^2$ matrix — a memory-efficient hierarchical representation of a data-sparse matrix whose row and column bases are nested across levels — can be constructed from scratch in time and memory proportional to the matrix dimension $N$, given only a black-box sketching routine (a fast matrix–vector product with random vectors) and a routine to evaluate matrix entries. Earlier sketching-based constructions worked top-down and consumed many random vectors, especially for three-dimensional problems, which made them slow and memory-hungry. The proposed algorithm works bottom-up: at each level of the cluster tree it strips off the contributions of already-known dense and coupling blocks from the sketched samples, applies a row interpolative decomposition to the residual samples to select skeleton rows, and propagates the reduced samples upward, adding fresh random vectors only when a local QR-based test says the node has not yet converged. On a GPU the same algorithm runs with batched kernels that handle variable-size blocks, and the paper reports up to $13\\times$ speedup over its own CPU version, up to $1000\\times$ speedup over an existing GPU top-down construction, and up to $660\\times$ speedup over an existing CPU sketching-based $\\mathcal{H}$ construction. If these claims hold, $\\mathcal{H}^2$ construction ceases to be the bottleneck it has been for large kernel matrices, integral equation operators, and low-rank updates inside sparse solvers.","feed_headline":"Adaptive GPU sketching builds H2 matrices 1000x faster","feed_subtitle":"A bottom-up algorithm cuts random-vector count from thousands to 256, unlocking larger dense problems on one GPU.","key_machinery":"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.","core_discovery":"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)$.","pith_inferences":["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."],"forward_implications":["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."],"supporting_citations":[{"why":"Supplies the HSS sketching and skeletonization machinery that Algorithm 1 extends from weak admissibility to strong admissibility and to adaptive sampling.","marker":"[29]"},{"why":"Earlier demonstration of randomized HSS construction inside multifrontal solvers; the application path the paper builds toward.","marker":"[7]"},{"why":"Top-down $\\mathcal{H}^2$ construction from matrix-vector products; the algorithmic baseline whose sample count grows with problem size.","marker":"[22]"},{"why":"GPU implementation of a top-down sketching-based $\\mathcal{H}^2$ construction used as the main GPU performance baseline.","marker":"[17]"},{"why":"Top-down $\\mathcal{H}$ construction with graph coloring used as the comparison baseline with $O(\\log N)$ sample count.","marker":"[23]"},{"why":"Fully black-box bottom-up recursive skeletonization; the class of algorithms this work brings to GPUs for the first time.","marker":"[24]"},{"why":"Supplies the marshaling and gathering primitives used to move data between flattened tree levels and batched kernels, with CPU and GPU backends.","marker":"[30]"},{"why":"Supplies batched dense linear algebra kernels used for the batched QR, gemm, and related operations.","marker":"[31]"},{"why":"Supplies batched factorization routines used for the non-uniform block-sparse row matrix multiplications.","marker":"[32]"}],"fun_headline_variants":["Bottom-up GPU sketching builds H2 matrices 1000x faster","First GPU bottom-up H2 construction: 1000x speedup","Adaptive sketching on GPU: H2 in O(N), 1000x faster","256 samples suffice: GPU H2 construction at 1000x","GPU H2 build with adaptive sketching: linear complexity"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Bottom-up GPU sketching builds H2 matrices 1000x faster","First GPU bottom-up H2 construction: 1000x speedup","Adaptive sketching on GPU: H2 in O(N), 1000x faster","256 samples suffice: GPU H2 construction at 1000x","GPU H2 build with adaptive sketching: linear complexity"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000828,"raw_usage":{"total_tokens":3666,"prompt_tokens":1041,"completion_tokens":2625,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":657,"completion_tokens_details":{"reasoning_tokens":2532}},"tokens_in":657,"tokens_out":2625,"duration_ms":19610,"temperature":1.0,"reasoning_tokens":2532,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T19:19:09.231873+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[{"cited_title":"A fast randomized algorithm for computing a hier- archically semiseparable representation of a matrix,","cited_arxiv_id":null,"evidence_quote":"Supplies the HSS sketching and skeletonization machinery that Algorithm 1 extends from weak admissibility to strong admissibility and to adaptive sampling."},{"cited_title":"A robust parallel preconditioner for indefinite systems using hierarchical matrices and randomized sampling,","cited_arxiv_id":null,"evidence_quote":"Earlier demonstration of randomized HSS construction inside multifrontal solvers; the application path the paper builds toward."},{"cited_title":"Fast construction of hierarchical matrix representation from matrix-vector multiplication,","cited_arxiv_id":null,"evidence_quote":"Top-down $\\mathcal{H}^2$ construction from matrix-vector products; the algorithmic baseline whose sample count grows with problem size."},{"cited_title":"H2Opus: a distributed-memory multi-GPU software package for non- local operators,","cited_arxiv_id":null,"evidence_quote":"GPU implementation of a top-down sketching-based $\\mathcal{H}^2$ construction used as the main GPU performance baseline."},{"cited_title":"Randomized Compression of Rank-Structured Matrices Accelerated with Graph Coloring","cited_arxiv_id":"2205.03406","evidence_quote":"Top-down $\\mathcal{H}$ construction with graph coloring used as the comparison baseline with $O(\\log N)$ sample count."},{"cited_title":"Randomized strong recursive skele- tonization: Simultaneous compression and factorization ofH-matrices in the black-box setting,","cited_arxiv_id":null,"evidence_quote":"Fully black-box bottom-up recursive skeletonization; the class of algorithms this work brings to GPUs for the first time."},{"cited_title":"Thrust: A productivity-oriented library for CUDA,","cited_arxiv_id":null,"evidence_quote":"Supplies the marshaling and gathering primitives used to move data between flattened tree levels and batched kernels, with CPU and GPU backends."},{"cited_title":"KBLAS: An optimized library for dense matrix-vector multiplication on GPU accelerators,","cited_arxiv_id":null,"evidence_quote":"Supplies batched dense linear algebra kernels used for the batched QR, gemm, and related operations."},{"cited_title":"Frame- work for batched and GPU-resident factorization algorithms to block Householder transformations,","cited_arxiv_id":null,"evidence_quote":"Supplies batched factorization routines used for the non-uniform block-sparse row matrix multiplications."}],"review_version":2}