{"id":"096ec3f0-ad9d-4cb5-a7f9-3b90aa142aa0","arxiv_id":"2509.11152","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":7,"one_line_summary":"An H2 strong-admissibility skeletonization solver with QR-based basis augmentation that factors and solves dense structured systems in linear time and memory, demonstrated on CPUs.","lead":"This paper presents a direct solver for dense linear systems with H2 hierarchical low-rank structure whose factorization and solve phases scale linearly in time and memory. It is aimed at making such solvers practical on modern many-core and GPU-style hardware for integral equations, Gaussian process kernels, and kernel methods.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Linear-complexity claim rests on an unproven bounded-rank premise that the paper itself flags as a 'sweet spot' in §3.1; if fill-in ranks grow with n, the O(n) bounds in §2.4 collapse.","rationale":"The reader's weakest assumption is the same as the load-bearing concern I find: the linear complexity proof is conditional on a bounded rank that is not established. I considered other plausible concerns—the absence of GPU measurements for a paper claiming fine-grained parallel architectures, the missing code artifact despite a release statement, and the lack of baseline comparisons. These are real scope/reproducibility limitations, but they do not threaten the correctness of the CPU solver's demonstrated scaling. The bounded-rank premise is more fundamental: if it fails, the central algorithmic complexity claim collapses even in the CPU-only setting, because the factorization and solve costs in §2.4 scale polynomially in k. The paper's own §3.1 caveat ('sweet spot', 'some intermediate recompression may be required') is an explicit admission that the demonstrated problems do not establish the general theorem. The paper deserves credit for honest reporting of this limitation, but the title and abstract state linear complexity without the qualification. The proposed test—tracking augmented ranks and time exponents at larger n—would determine whether the concern is real for the presented families or only a generality gap. If the test shows rank growth, the central claim should be downgraded; if not, the conditional verdict can stand or be upgraded. Since the reader already arrived at CONDITIONAL with high confidence, no verdict change is needed.","tokens_in":15159,"tokens_out":14661,"duration_ms":178247,"concrete_test":"Instrument the released CPU solver to log, for the 3D Helmholtz and 3D covariance cases of §3.1, the maximum augmented basis rank k_aug_max at each level as n is increased from 2^14 to 2^20 and, if feasible, to 2^22. Fit log(k_aug_max) vs log(n) and also fit the log of factorization time vs log(n). If k_aug_max grows with n (slope > 0) or the time exponent exceeds 1, the bounded-rank premise fails for these families and the observed linearity is a pre-asymptotic artifact; if k_aug_max stays flat, the concern is mitigated for the demonstrated problems.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim of a linear-complexity H2 direct solver depends on the premise, stated in §2.4, that the rank k—including the ranks of augmented bases for fill-in—is bounded by a constant independent of n. This premise is not proved. Section 2.1 (step 4) asserts only that fill-in blocks are 'admissible'; admissibility guarantees low-rank approximability, not a uniform bound on the approximating rank as Schur complements accumulate. Section 3.1 concedes that the four test families are a 'sweet spot' in which 'block ranks do not grow substantially during the factorization,' and that 'in more general cases... some intermediate recompression may be required.' Thus the theorem-like headline 'Linear Complexity' is conditional: the measured linear scaling up to n=2^20 cannot distinguish a genuinely linear algorithm from one with slowly growing rank (e.g., k=O(log n) or O(n^ε)), whose complexity would be superlinear. Since recompression is not part of the presented algorithm, the claim as stated overreaches the evidence.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents factorization and solve phases for a strong-admissibility H2 direct solver based on recursive skeletonization, with a focus on batched, fine-grained parallel execution. The algorithmic contributions are a multi-level graph coloring to identify independent skeletonization steps, QR-based basis augmentation to accommodate Schur-complement fill-in, prefix-sum memory management, and a hierarchical forward/backward solve. The complexity analysis in §2.4 claims O(C_sp^2 k^3 n) factorization and O(C_sp k^2 n) solve, linear under a bounded-rank assumption. Experiments on four dense matrix families (2D/3D covariance, 2D Laplace integral equation, 3D Helmholtz integral equation, low-rank update) up to n=2^20 show near-linear time and memory scaling and backward errors near the requested tolerance. A CPU implementation is measured; a GPU implementation is stated as in progress.","tokens_in":15465,"tokens_out":4798,"duration_ms":59998,"significance":"If the linear-complexity claim holds under the stated assumptions, the solver would be a valuable tool for dense and data-sparse systems, particularly because it is blackbox with respect to the origin of the matrix and is designed for batched execution. The paper gives unusually detailed algorithmic steps and includes an experimental backward-error study, and the authors state an intention to release code. However, the central complexity theorem is explicitly conditional on a uniform bound on numerical rank k, including ranks of fill-in bases, and this premise is not proved. The experiments are deliberately in a 'sweet spot' where ranks do not grow substantially. The unqualified 'linear complexity' claim in the title and abstract therefore overreaches the evidence, and additional analysis or careful qualification is needed before the central claim can be accepted.","major_comments":[{"comment":"The complexity bound is stated in §2.4 as 'if we assume a bound k on the rank.' This is precisely the load-bearing assumption. The manuscript does not prove that the augmented bases or the ranks of fill-in blocks generated during Schur complement updates remain bounded by a constant independent of n. Section 3.1 concedes that the four test problems 'occupy a sweet spot' and that 'in more general cases ... some intermediate recompression may be required.' Since recompression is not part of the presented algorithm, the advertised 'linear complexity direct solver' is conditional rather than demonstrated. The experiments up to n=2^20 cannot distinguish a truly linear algorithm from one with slowly growing rank, e.g., k=O(log n) or O(n^ε), whose complexity is superlinear. The authors should either provide a theorem with precise conditions under which the rank bound holds, or explicitly qualif","section":"§2.4 and §3.1"},{"comment":"The statement 'By construction, the fill-in blocks are always admissible' is used to justify that fill-in is low-rank, but admissibility is a geometric separation condition, not a rank bound. The augmented basis computed in Algorithm 2 via a truncated SVD of (I - V_i V_i^T) F_{i*} controls the one-step approximation error for a single cluster; it does not control the ranks of blocks generated by later partial LU updates, which can accumulate contributions from multiple clusters. The argument needs an induction establishing a uniform bound on the ranks of all fill-in blocks across all levels. This missing step is load-bearing for the O(n) factorization and solve bounds.","section":"§2.1, step 4"},{"comment":"The evidence for linear scaling is visual: O(n) reference lines are overlaid on log-log plots, but no slopes, confidence intervals, or numbers of runs are reported. With the tested range up to n=2^20, the plots cannot rule out log-linear or mildly superlinear behavior. I recommend reporting factorization and solve times in tabular form, fitting slopes with confidence intervals, plotting the maximum rank or rank distribution as a function of n for each family, and comparing against a baseline solver (e.g., a standard H2 or HODLR solver). Also state whether the timings are single runs or averaged, and include error bars if repeats were performed.","section":"§3 and Figs. 13, 16"},{"comment":"The thread-scaling results are not uniformly consistent with the abstract's claim of 'parallel scaling up to 16 threads.' For the 2D factorization, the time at 16 threads (7.881 s) is worse than at 8 threads (7.099 s), and the large QR batched operation also degrades between 8 and 16 threads (179.65 ms to 198.18 ms). The text acknowledges 'seemingly poor scaling' but the abstract's wording is stronger than the data. Please either qualify the scaling claim or investigate the source of the degradation (e.g., memory bandwidth, false sharing, OpenMP overhead) and report whether the same behavior would be expected on a GPU.","section":"§3.3, Table 3"}],"minor_comments":[{"comment":"Typo: 'motivated by the needed to reduce rank growth' should be 'motivated by the need to reduce rank growth.'","section":"§1.2"},{"comment":"Minor formatting: 'correlations lengthslset to' should read 'correlation lengths l set to'.","section":"§3.1"},{"comment":"Line 6, 'F=fillIn(A,F,l)', is not defined in the pseudocode or text. Clarify how the fill-in matrix is extracted and how it interacts with levels.","section":"Algorithm 1"},{"comment":"Line 3 uses the notation V_i for the cluster basis and F_{i*} for a block row of the fill-in matrix, but the dimensions are not stated. Add a sentence explaining the exact matrix view used in the SVD, especially since the fill-in blocks may be rectangular.","section":"Algorithm 2"},{"comment":"The paper says 'we are releasing code for the solver in open-source,' but no repository URL or availability statement is given. Please provide a link or a footnote with the intended release location.","section":"Conclusion"},{"comment":"Figures 13 and 16 should state the number of runs and whether the plotted points are medians, means, or single measurements. Also clarify what the O(n) reference lines mean (e.g., arbitrary constant offset).","section":"Figure captions"}],"recommendation":"major_revision","confidential_remarks":"The paper is a solid algorithmic contribution, but the headline claim overreaches the theoretical support. In addition to the bounded-rank issue, the experimental section needs quantitative evidence and an honest treatment of thread scaling. The novelty relative to Ma et al. and Ma/Yokota should be sharpened in revision, particularly the claimed benefit of QR-based augmentation and the GPU-oriented batched formulation, since no GPU results are actually presented. A major revision with a qualified central claim and additional experiments is appropriate."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Colleague,\n\nRead this if you work on hierarchical direct solvers or GPU batched linear algebra. It's a careful CPU implementation of a strong-admissibility H2 skeletonization solver. The genuinely new pieces are: QR-based basis augmentation for fill-in (avoids the condition-number squaring of Gram approaches), multi-level graph coloring for conflict-free batching, prefix-sum memory management, and a parallel solve phase. The experiments cover four problem families up to n=2^20, backward error is measured and hits the requested thresholds, and the phase-by-phase timing breakdown is honest about memory-bandwidth limits. The complexity analysis in §2.4 is standard and clearly stated.\n\nThe main soft spot is the linear-complexity headline. The O(n) bounds assume a constant rank bound k, including fill-in ranks. That premise is asserted, not proved, and §3.1 concedes the test problems are a 'sweet spot' and that general blackbox use may require recompression. So the title overreaches, though this is a common assumption in the H2 literature and not fatal to the paper's practical value.\n\nOther soft spots: GPU-centric framing with no GPU data—the results are on 16 CPU threads. No baseline against existing solvers like Ma-Yokota or FMM-LU, so prefactors aren't positioned. No error bars or repeated runs, though the scaling trends are clear. The text promises an open-source release but no artifact link appears in the manuscript.\n\nThe citation pattern is fair; the debt to Ma-Jiao and Ma-Yokota is explicit. The reader's circularity concern is correctly low: this is measurement, not parameter-fitting.\n\nVerdict: worth a serious referee. A reviewer should push for a precise statement of the rank-growth assumption, GPU results or a clear explanation of their absence, a baseline comparison, and the promised code. I'd accept it for review and expect major revision.","headline":"Solid engineering paper with honest experiments, but the 'Linear Complexity' title outruns the bounded-rank assumption the authors themselves flag as a sweet spot.","tokens_in":15931,"tokens_out":2430,"would_cite":true,"duration_ms":26685,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65F05","65F30","65Y05"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper presents a blackbox direct solver for strong-admissibility H2 matrices that factors and solves dense structured systems in linear time and linear memory, demonstrated on systems up to one million unknowns.","keywords":["H2 matrices","direct linear solver","strong admissibility","recursive skeletonization","linear complexity","hierarchical low-rank matrices","parallel computation","backward error"],"falsifier":"Construct a strongly admissible H2 matrix whose off-diagonal blocks have rank that grows with the matrix size—for example, a Helmholtz kernel with wavenumber proportional to the mesh resolution, or a covariance matrix with a short correlation length—factor it at increasing sizes, and check whether the time per unknown and memory per unknown stay flat. If either grows with n, the linear-complexity claim is falsified for that class of problems.","tokens_in":15073,"feed_emoji":"🧮","tokens_out":4689,"duration_ms":50129,"temperature":0.7,"pith_summary":"The paper aims to establish that strong-admissibility H2 matrices—dense structured matrices whose far-field interactions are compressed into nested low-rank bases—can be factored and solved by a direct, non-iterative solver in time and memory that scale linearly with the number of unknowns. On four families of dense systems, including covariance matrices, Laplace and Helmholtz integral operators, and a global low-rank update, it reports factorization and solve timings and memory usage that track the linear reference lines up to 2^20 unknowns, with backward error controlled by the chosen truncation thresholds. The practical payoff would be a robust direct solver for large dense systems from integral equations, Gaussian processes, and related applications, designed for the batched, fine-grained parallelism of GPUs and many-core CPUs.","feed_headline":"Linear-time direct solver reaches one-million-unknown dense systems","feed_subtitle":"Hierarchically low-rank dense systems from covariance models and integral equations now factor in linear time and memory.","key_machinery":"The key machinery is the strong-admissibility H2 format, where admissible blocks share nested basis matrices and small coupling matrices, combined with a fill-in matrix F that is processed together with the admissible part H. On each cluster, the algorithm (i) forms an augmented orthogonal basis spanning the existing cluster basis and the dominant directions of the fill-in block row, using a truncated SVD of a QR factor; (ii) applies the full orthogonal projection to the admissible, inadmissible, and fill-in parts; and (iii) performs a partial LU elimination on the inadmissible diagonal block, moving the resulting admissible fill-in into F. Multi-level graph coloring partitions clusters into","core_discovery":"The central claim is that recursive skeletonization, specialized to the strong-admissibility H2 format, yields a factorization expressed as a product of block unit-triangular factors and orthogonal projections, and that both factorization and solve execute in strictly linear time and memory. The algorithm augments the nested cluster bases on the fly with a QR-based basis for the fill-in produced by Schur complement updates, so that fill-in blocks remain admissible and are skeletonized together with the original admissible blocks. The measured complexity is O(Csp^2 k^3 n) for factorization and O(Csp k^2 n) for the solve, which is linear in n whenever the numerical rank k and the sparsity cons","pith_inferences":["Inference: The paper's linear-complexity claim rests on an unproven bounded-rank assumption, and its own test problems are described as a 'sweet spot.' For problems with genuinely rank-growing blocks—such as high-frequency oscillatory kernels—the algorithm as presented may need intermediate recompression before it is truly blackbox for all H2 matrices.","Inference: The QR-then-SVD basis augmentation avoids the condition-number squaring of Gram-matrix approaches, so the method is structurally well suited to lower-precision execution; the paper discusses lower-precision alternatives but does not test them, making single- or mixed-precision runs a natural next experiment.","Inference: The practical constants in the O(Csp^2 k^3 n) estimate are strongly affected by problem dimension and by choices such as leaf size and admissibility parameter η; users may need to tune these per application, especially for 3D problems where observed sparsity constants and ranks are substantially larger."],"forward_implications":["If the reported scaling holds, direct solves of dense structured systems with one million unknowns are practical on a single many-core CPU node, with memory usage growing only linearly in n.","The blackbox nature—only the H2 matrix and right-hand side are needed, without geometric or analytic information—means the solver can be dropped into existing computational workflows for Gaussian processes, integral equations, and kernel methods.","The solve phase shares the hierarchical matrix-vector structure, so forward and backward substitutions for multiple right-hand sides inherit the same linear-time per-solve cost.","Backward errors stay at the level of the specified truncation thresholds (relative errors in the 1e-6 to 1e-7 range in the reported tests), supporting use as a numerical solver rather than merely a preconditioner.","Because all operations are marshaled into batches with preplanned memory allocation, the algorithm is structured for GPU execution; the paper reports a CPU implementation and states that a GPU implementation is in progress."],"fun_headline_variants":["H2 solver factors million-size dense systems in linear time","Linear-scaling direct solver for dense matrices up to 1M unknowns","Direct H2 solver: linear time and memory to a million unknowns","Skeletonized H2 factorization achieves linear cost for dense solves","One-million dense system solved in linear time with H2"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The entire linear-time argument assumes that the numerical rank of every admissible block and every fill-in block stays below a constant that is independent of the problem size n; if ranks grow with n (for example, with increasing wave number or shrinking correlation length), the O(n) complexity no longer follows and the factorization would require recompression.","fun_headline_variants_meta":{"raw":{"variants":["H2 solver factors million-size dense systems in linear time","Linear-scaling direct solver for dense matrices up to 1M unknowns","Direct H2 solver: linear time and memory to a million unknowns","Skeletonized H2 factorization achieves linear cost for dense solves","One-million dense system solved in linear time with H2"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000598,"raw_usage":{"total_tokens":2621,"prompt_tokens":723,"completion_tokens":1898,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":467,"completion_tokens_details":{"reasoning_tokens":1811}},"tokens_in":467,"tokens_out":1898,"duration_ms":16975,"temperature":1.0,"reasoning_tokens":1811,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-04T16:59:52.810869+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Construct a strongly admissible H2 matrix whose off-diagonal blocks have rank that grows with the matrix size—for example, a Helmholtz kernel with wavenumber proportional to the mesh resolution, or a covariance matrix with a short correlation length—factor it at increasing sizes, and check whether the time per unknown and memory per unknown stay flat. If either grows with n, the linear-complexity claim is falsified for that class of problems.","supporting_citations":[],"review_version":1}