REVIEW 4 major objections 6 minor 39 references
Hierarchical Recursive Precision for Accelerating Symmetric Linear Solves on MXUs
T0 review · 4 major / 6 minor · reviewed 2026-08-03 · deepseek-v4-flash
Pith's one-line read A fully nested recursive Cholesky factorization can assign low-precision FP16 arithmetic to large off-diagonal blocks while keeping diagonal blocks in higher precision, yielding multi-fold speedups on AI-specialized GPUs without sacrificing
desk verdict Plausible speedups and a useful portable design, but the stability claim is narrower than advertised and the evidence needs an artifact release. 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 recursive decomposition tree: at each level the matrix is split into a diagonal block and a trailing submatrix, and all three core operations — diagonal factorization (POTRF), triangular solve (TRSM), and symmetric rank-k update (SYRK) — are applied recursively rather than in standard blocked form. This tree is the data structure that carries the mixed-precision hierarchy: precision is assigned per recursion level, with large off-diagonal blocks handled by FP16 GEMM and diagonal blocks by FP32/FP64. A per-block quantization step rescales blocks into FP16's representable range before low-precision operations and rescales them back afterward, which the paper argues is necessary and cheap.
What would settle it
Measure the factorization error ||A - LL^T||/||A|| for the deepest FP16 configuration on a dense SPD matrix with condition number 10^6 and n=65,536; the paper predicts roughly 5–6 digits, while pure FP16 gives below 4. If the measured error matches pure FP16 (below 4 digits), the claimed accuracy benefit is falsified. Alternatively, a single run on a matrix with dynamic range above 1e11 (e.g., raefsky4) producing NaN would show the static scaling assumption fails.
Extended reading notes
Core claim
The central discovery is that the numerical sensitivity of Cholesky factorization is concentrated near the diagonal for diagonally dominant SPD matrices, so the abundant off-diagonal GEMM work can be executed in FP16 without destroying stability, provided diagonal blocks are computed in higher precision and each low-precision block is rescaled to FP16's range. The paper implements this as a recursive tree that splits the matrix at each level into two recursive factorizations, one triangular solve, and one symmetric rank-k update, assigning lower precision to blocks further from the diagonal. The claimed payoff is quantitative: up to 5.32x speedup over the vendor FP64 baseline on H200 with 5–
Load-bearing premise
For diagonally dominant SPD matrices, the off-diagonal blocks are insensitive enough that static FP16 scaling does not destroy the factorization, and this holds as recursion depth increases.
Editorial extensions
If this is right
- Large diagonally dominant SPD systems can be solved several times faster than double-precision vendor libraries on existing tensor-core hardware.
- Recursive SYRK and fully nested recursive mixed-precision Cholesky are feasible on both NVIDIA and AMD GPUs from a single portable implementation.
- Accuracy can be traded continuously with speed by choosing recursion depth and precision layering (e.g., a configuration with a top FP64 layer yields roughly 9 digits at 1.21x speedup, while deeper FP16 layers yield 5–6 digits at 5.32x speedup).
- The approach is least effective for small matrices (n < 16,000) and fails for matrices with extreme dynamic range, where static blockwise quantization causes underflow and NaN propagation.
- Peak utilization is not the right objective: layered mixed precision can trade a modest loss of raw throughput for dramatically improved accuracy while still delivering multi-fold speedup over high-precision solvers.
Reading between the lines
- The same tree-structured precision hierarchy could extend to LDL^T factorizations and indefinite systems, where diagonal pivoting complicates precision assignment but the off-diagonal insensitivity argument may partially carry over.
- Because accuracy is measured as factorization error rather than solve residual, the 5–6 digit claim may be conservative for well-conditioned right-hand sides; combining this solver with iterative refinement could recover full precision at modest extra cost, making the speedup usable in production.
- The static blockwise scaling is the clearest extension point: adaptive or per-block dynamic scaling based on block norms or exponent ranges might extend the method to matrices like raefsky4 without losing the performance advantage.
- Performance portability across AMD and NVIDIA depends on the availability of mixed-precision GEMM; once AMD's mixed-precision GEMM matures, the AMD speedups may approach the NVIDIA numbers.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a nested recursive Cholesky factorization for symmetric positive-definite systems, in which not only POTRF but also TRSM and SYRK are recursively subdivided. A layered mixed-precision scheme assigns FP16 to large off-diagonal blocks and higher precision (FP32/FP64) to diagonal blocks, with per-block quantization to mitigate FP16 dynamic-range issues. The implementation is in Julia and is evaluated on NVIDIA H200 and AMD MI300X. The central quantitative claim is up to a 5.32x speedup over cuSOLVER FP64 at n=65,536 while retaining 5-6 digits of factorization accuracy, about 100x better than pure FP16. The paper also claims the first recursive GPU SYRK and cross-vendor portability.
Significance. If the numerical and performance claims hold, the work is a useful contribution: it combines recursive blocking with tensor-core/matrix-core low-precision arithmetic in a portable way, and it addresses a real bottleneck in symmetric direct solvers. The Julia multiple-dispatch design is a credible route to hardware-agnostic mixed-precision linear algebra. The authors are transparent about several limitations (small matrices, dynamic-range sensitivity, AMD missing mixed-precision GEMM), which increases confidence in the reported experiments. However, the central numerical-stability premise is not proven for the stated target class, and the performance results are not accompanied by code, data, or error bars, so reproduction is not currently possible.
major comments (4)
- [§3.2, §4.1–4.2, §4.5] The numerical-stability premise is load-bearing but unproven. §3.2 asserts that for diagonally dominant matrices 'the most numerically sensitive operations' are near the diagonal, but the benchmarks use A = rand(n) + nI, which has condition number near 1, and one moderately conditioned real matrix (bodyy5). A strictly diagonally dominant SPD matrix with O(1) entries can have condition growing with n (e.g., the 2D Laplacian); FP16 off-diagonal rounding (~1e-3 relative) may then be amplified through the Schur complement, and the claimed 5–6 digits of accuracy are not established. The failure reported for raefsky4 in §4.5 is attributed to dynamic range, but dynamic range is not the only possible failure mode. Without an error bound or experiments on diagonally dominant matrices with growing condition number, the headline accuracy claim is unsupported for the stated target class.
- [Figures 4–12; §4.1] The performance claims are not reproducible from the manuscript: no code, no data, no per-run statistics, and no error bars are provided. The vendor baselines (cuSOLVER FP64/FP32, rocSOLVER FP64) and any tuning parameters for the recursive kernels or leaf block sizes are not described in enough detail for an independent implementation. For a systems/performance paper, this is a major gap: the central speedup figures cannot currently be verified.
- [§4.6, Figure 12] The cross-platform comparison is not apples-to-apples. The text states that mixed-precision GEMM (GemmEx) was included for NVIDIA but not for AMD because of availability in the Julia ecosystem. Thus the lower MI300X speedups and the portability claim in Figure 12 conflate algorithmic portability with the absence of an available vendor path on one platform. The claim of 'similar gains' on AMD is therefore weaker than the NVIDIA results.
- [Abstract and §4.2] The headline performance number is reported inconsistently. The abstract says '5.07x speedup relative to the diagonal-precision vendor baseline' (presumably 'double-precision'), while §4.2 and the introduction report 5.32x versus cuSOLVER FP64. The deepest mixed-precision configuration is also written with different numbers of FP16 layers in the same section: [F16,F16,F16,F16,F16,F32], [F16,F16,F16,F16,F16,F16,F32], and [F16,F16,F16,F16,F16,F16,F16,F32] are all mentioned. This ambiguity needs to be resolved because the 5.32x claim is the central result.
minor comments (6)
- [Abstract] 'diagonal-precision' appears to be a typo for 'double-precision'; please correct.
- [§3.3] The scaling formula defines alpha but does not explicitly state that B_alg is set to B/alpha before the low-precision operation and that the output is multiplied by alpha; make this two-step explicit.
- [§2.1] The claim of 'first recursive GPU-based SYRK' would be easier to evaluate if the related-work search were explicit; currently no comparison to any prior recursive SYRK implementation is provided.
- [§4.5] The text refers to Figure 11 when discussing small-matrix overhead, but Figure 11 is a scaling plot; consider referencing a more appropriate figure or clarifying the intended comparison.
- [Figure 8] The y-axis label '-log10 relative error' can be misread as 'digits of accuracy'. Define the relationship explicitly.
- [Algorithm 1] The partition uses A21 in the split but the algorithm description says lower Cholesky; check that the transpose notation is consistent with a lower-triangular factor.
Circularity Check
No significant circularity; empirical benchmarks and accuracy comparisons are self-contained.
full rationale
The paper's central claims are quantitative benchmark results: mixed-precision recursive Cholesky configurations achieve up to 5.32x speedup over cuSOLVER FP64 with 5-6 digits of factorization accuracy. These numbers are measured against external vendor baselines (cuSOLVER/rocSOLVER) and an independent FP64 factorization residual, not derived from the method's own assumptions. There is no fitted parameter that is subsequently renamed as a prediction, and no equation in Sections 3.1-3.3 defines the claimed accuracy or speedup in terms of the inputs. The quantization/quantization scheme is a deterministic scaling operation, not a model fit. The numerical-sensitivity assumption for diagonally dominant matrices is stated qualitatively in Section 3.2 and is not proven, but an unproven assumption is a correctness/stability risk, not circularity. The paper explicitly documents the failure on Simon/raefsky4, so it does not hide its limitation. The only self-citations ([16], [17]) refer to the authors' earlier recursive TRSM work and a poster; these are used as background and as evidence of portable Julia GPU kernels, but the current paper's own benchmarks and external baselines carry the central argument. No uniqueness theorem or imported ansatz is invoked to force the mixed-precision hierarchy. Therefore no circular step can be exhibited with a specific reduction, and the appropriate finding is no significant circularity, with a minor non-load-bearing self-citation noted.
Assumptions & free parameters
free parameters (3)
- leaf block size b =
not reported
- precision hierarchy configuration =
various, e.g. [F16,F16,F16,F16,F16,F16,F32]
- diagonal shift n for synthetic matrices =
n (matrix size)
assumptions (4)
- standard math Recursive block decomposition of an SPD matrix into diagonal and off-diagonal blocks yields the same Cholesky factors as the standard algorithm.
- domain assumption Off-diagonal blocks are numerically less sensitive than diagonal blocks for diagonally dominant SPD matrices, so FP16 round-off there is acceptable.
- domain assumption Static per-block quantization scaling alpha keeps values in FP16 range without materially changing the result for matrices of limited dynamic range.
- domain assumption Vendor BLAS/LAPACK kernels (POTRF, TRSM, SYRK, GEMM) used for base cases meet the expected accuracy and performance on both H200 and MI300X.
Cite this review
Pith. "Pith review of Hierarchical Recursive Precision for Accelerating Symmetric Linear Solves on MXUs." pith.science (2026). https://pith.science/paper/Y2JTYCTI
@misc{pith2026260108082,
author = {Pith},
title = {Pith review of: Hierarchical Recursive Precision for Accelerating Symmetric Linear Solves on MXUs},
year = {2026},
howpublished = {\url{https://pith.science/paper/Y2JTYCTI}},
note = {Machine review of arXiv:2601.08082}
}
read the original abstract
Symmetric positive-definite system solvers based on Cholesky factorization are fundamental to many scientific applications, such as climate modeling. We present a portable, nested recursive mixed-precision solver designed for Matrix Processing Units (MXUs), including NVIDIA Tensor Cores (H200) and AMD Matrix Cores (MI300X), that assigns low-precision FP16 arithmetic to large off-diagonal blocks, while preserving high precision on diagonal blocks to ensure numerical stability. The solver is implemented in Julia, providing a high-level, hardware-agnostic interface. We demonstrate up to a 5.07x speedup relative to the diagonal-precision vendor baseline, with 100x better accuracy than pure half precision on H200, providing higher accuracy than low-precision at higher speed than high-precision. Positive performance trends are also observed on MI300X, demonstrating broad applicability across GPUs.
Figures
Figures from the paper (5 more)
Reference graph
Works this paper leans on
-
[1]
Agullo, E., Beaumont, O., Eyraud-Dubois, L., et al.: Bridging the gap between performance and bounds of Cholesky factorization on het- erogeneous platforms. In: Proc. IPDPS Workshops. pp. 34–45 (2015). https://doi.org/10.1109/IPDPSW.2015.35
-
[2]
In: Proc
Akbudak, K., Ltaief, H., Mikhalev, A., Keyes, D.: Tile low rank Cholesky factor- ization for climate/weather modeling applications on manycore architectures. In: Proc. ISC High Performance. pp. 22–40. Springer (2017)
2017
-
[3]
In: Proc
Alomairy, R., Abdulah, S., Cao, Q., Genton, M.G., Keyes, D.E., Ltaief, H.: Sus- tainably modeling a sustainable future climate. In: Proc. HPEC. pp. 1–8. IEEE (2025)
2025
-
[4]
In: Proc
Alomairy, R., Cao, Q., Ltaief, H., Keyes, D., Edelman, A.: Scalable Hamming distance computation using accelerated matrix transformations. In: Proc. ISC High Performance. pp. 1–13 (2025)
2025
-
[5]
In: Proc
Alomairy, R., Tome, F., Samaroo, J., Edelman, A.: Dynamic task scheduling with data dependency awareness using Julia. In: Proc. HPEC. pp. 1–7. IEEE (2024)
2024
-
[6]
Alomairy, R.M.: High-performance scientific applications using mixed precision and low-rank approximation powered by task-based runtime systems. Ph.D. thesis, KAUST (2022)
2022
-
[7]
Amestoy, P.R., et al.: Mixed precision low-rank approximations and their applica- tion to block low-rank LU factorization. IMA J. Numer. Anal.43(4), 2198–2227 (2023) 14 V. Carrica et al
2023
-
[8]
ACM Trans
Andersen, B.S., Wa´ sniewski, J., Gustavson, F.G.: A recursive formulation of Cholesky factorization of a matrix in packed storage. ACM Trans. Math. Softw. 27(2), 214–244 (2001)
2001
Show all 39 references
-
[9]
SIAM (1999)
Anderson, E., et al.: LAPACK Users’ Guide. SIAM (1999)
1999
-
[10]
Computer Physics Communications180(12), 2526–2533 (2009)
Baboulin, M., et al.: Accelerating scientific computations with mixed precision algorithms. Computer Physics Communications180(12), 2526–2533 (2009)
2009
-
[11]
In: Proc
Ballard, G., Demmel, J., Holtz, O., Schwartz, O.: Communication-optimal parallel and sequential Cholesky decomposition. In: Proc. SPAA ’09. pp. 245–252 (2009). https://doi.org/10.1145/1583991.1584054
2009
-
[12]
SIAM Review59(1), 65–98 (2017)
Bezanson, J., Edelman, A., Karpinski, S., Shah, V.B.: Julia: A fresh approach to numerical computing. SIAM Review59(1), 65–98 (2017)
2017
-
[13]
Computer Science – Research and Devel- opment29(2), 139–147 (2014)
Bosilca, G., Ltaief, H., Dongarra, J.: Power profiling of Cholesky and QR factor- izations on distributed memory systems. Computer Science – Research and Devel- opment29(2), 139–147 (2014)
2014
-
[14]
Buttari, A., et al.: Mixed precision iterative refinement techniques for the solution of dense linear systems. Int. J. High Perform. Comput. Appl.21(4), 457–466 (2007)
2007
-
[15]
In: Proc
Cao, Q., et al.: Reducing data motion and energy consumption of geospatial model- ing applications using automated precision conversion. In: Proc. IEEE CLUSTER (2023)
2023
-
[16]
Carrica, V., Alomairy, R., Ringoot, E., Edelman, A.: Accelerating linear solve with mixed precision nested recursive subdivision on AI hardware (2025), sC25 Poster
2025
-
[17]
arXiv preprint arXiv:2504.13821 (2025)
Carrica, V., Onyango, M., Alomairy, R., Ringoot, E., Schloss, J., Edelman, A.: Toward portable GPU performance: Julia recursive implementation of TRMM and TRSM. arXiv preprint arXiv:2504.13821 (2025)
2025 arXiv
-
[18]
Carson, E., Higham, N.J.: A new analysis of iterative refinement and its application to accurate solution of ill-conditioned sparse linear systems. SIAM J. Sci. Comput. 39(6), A2834–A2856 (2017)
2017
-
[19]
Carson, E., Higham, N.J.: Accelerating the solution of linear systems by iterative refinement in three precisions. SIAM J. Sci. Comput.40(2), A817–A847 (2018)
2018
-
[20]
Concurrency and Computation: Practice and Experience29(15), e4187 (2017)
Charara, A., Keyes, D., Ltaief, H.: A framework for dense triangular matrix kernels on various manycore architectures. Concurrency and Computation: Practice and Experience29(15), e4187 (2017)
2017
-
[21]
In: Proc
Chen, C., Martinsson, P.G.: Solving linear systems on a GPU with hierarchically off-diagonal low-rank approximations. In: Proc. SC22. pp. 1–15. IEEE (2022)
2022
-
[22]
Eliahu, D., Spillinger, O., Fox, A., Demmel, J.: Frpa: A framework for recursive parallel algorithms. Tech. rep., EECS, University of California at Berkeley (2015)
2015
-
[23]
Concurrency and Computation: Practice and Experience35(25), e7811 (2023)
Faverge, M., et al.: Programming heterogeneous architectures using hierarchical tasks. Concurrency and Computation: Practice and Experience35(25), e7811 (2023)
2023
-
[24]
Furmento, N., et al.: Optimizing parallel heterogeneous system efficiency: Dy- namic task graph adaptation with recursive tasks. J. Parallel Distrib. Comput. 205, 105157 (2025)
2025
-
[25]
Ad- vances in Neural Information Processing Systems31(2018)
Gardner, J., Pleiss, G., Weinberger, K.Q., Bindel, D., Wilson, A.G.: GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. Ad- vances in Neural Information Processing Systems31(2018)
2018
-
[26]
In: Proc
Giordano, M., Kl¨ ower, M., Churavy, V.: Productivity meets performance: Julia on A64FX. In: Proc. IEEE CLUSTER. pp. 549–555 (2022)
2022
-
[27]
Numerische Mathematik112(4), 565–600 (2009) Hierarchical Recursive Precision for POTRF on MXUs 15
Grasedyck, L., Kriemann, R., Le Borne, S.: Domain decomposition based H-LU preconditioning. Numerische Mathematik112(4), 565–600 (2009) Hierarchical Recursive Precision for POTRF on MXUs 15
2009
-
[28]
In: Proc
Haidar, A., Tomov, S., Dongarra, J., Higham, N.J.: Harnessing GPU tensor cores for fast FP16 arithmetic to speed up mixed-precision iterative refinement solvers. In: Proc. SC18. pp. 603–613. IEEE (2018)
2018
-
[29]
Acta Numerica31, 347–414 (2022)
Higham, N.J., Mary, T.: Mixed precision algorithms in numerical linear algebra. Acta Numerica31, 347–414 (2022)
2022
-
[30]
In: Proc
Lang, J., Guo, Z., Huang, S.: A comprehensive study on quantization techniques for large language models. In: Proc. ICAIRC. pp. 224–231. IEEE (2024)
2024
-
[31]
In: Proc
Ltaief, H., Alomairy, R., Cao, Q., et al.: Toward capturing genetic epistasis from multivariate genome-wide association studies using mixed-precision kernel ridge regression. In: Proc. SC24. pp. 1–12. IEEE (2024)
2024
-
[32]
In: Proc
Ltaief, H., Tomov, S., Nath, R., Du, P., Dongarra, J.: A scalable high performant Cholesky factorization for multicore with GPU accelerators. In: Proc. VECPAR
-
[33]
Aerospace10(9), 792 (2023)
Piscaglia, F., Ghioldi, F.: GPU acceleration of CFD simulations in OpenFOAM. Aerospace10(9), 792 (2023)
2023
-
[34]
In: ISC High Performance 2025 Research Paper Proceedings (40th International Conference)
Ren, J., Ltaief, H., Abdulah, S., Keyes, D.E.: Accelerating mixed-precision out-of- core cholesky factorization with static task scheduling. In: ISC High Performance 2025 Research Paper Proceedings (40th International Conference). pp. 1–12. Prom- eteus GmbH (2025)
2025
-
[35]
arXiv preprint arXiv:2510.12705 (2025)
Ringoot, E., Alomairy, R., Edelman, A.: A GPU-resident memory-aware al- gorithm for accelerating bidiagonalization of banded matrices. arXiv preprint arXiv:2510.12705 (2025)
2025 arXiv
-
[36]
Advanced Theory and Simulations8(3), 2400629 (2025)
Tomanovi´ c, I., et al.: CFD code parallelization on GPU and the code portability. Advanced Theory and Simulations8(3), 2400629 (2025)
2025
-
[37]
Com- puting in Science & Engineering11(6), 56–63 (2009)
Van Zee, F.G., et al.: The libflame library for dense matrix computations. Com- puting in Science & Engineering11(6), 56–63 (2009)
2009
-
[38]
In: Proc
Zhang, Q., Alomairy, R., Wang, D., Gu, Z., Cao, Q.: Leveraging hardware-aware computation in mixed-precision matrix multiply: A tile-centric approach. In: Proc. W AMTA. pp. 174–185. Springer (2025)
2025
-
[2010]
pp. 93–101. Springer (2011)
2011
Reviewed August 3, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.