This article presents two new multilevel hierarchical matrices to perform fast kernel matrix-vector products in $d$ dimensions $(d>1)$, namely $\mathcal{H}_{hw}^{\prime}$, $\mathcal{H}_{hw}^{\prime \prime}$. We introduce a novel hybrid weak admissibility condition to construct these hierarchical matrices. More precisely, at coarser levels of the tree (i.e., larger blocks), we use a weak admissibility condition, wherein \textbf{far-field} and \textbf{vertex-sharing} clusters are considered admissible, while at the finer levels (i.e., smaller blocks), the \textbf{adjacent} clusters are treated as admissible. The admissible blocks at coarser levels of $\mathcal{H}_{hw}^{\prime}$ and $\mathcal{H}_{hw}^{\prime \prime}$ are represented using semi-nested and fully nested bases, respectively. In both hierarchical matrices, the finer levels at which adjacent clusters are treated as admissible are constructed in a non-nested fashion. Thus, they can be viewed as a hybrid of $\mathcal{H}$ and $\mathcal{H}^2$ matrices. We use adaptive cross approximation (ACA) and nested cross approximation (NCA) to construct non-nested and nested bases, respectively. The fast matrix-vector products are performed using the proposed hierarchical matrix formats. To evaluate the performance and usefulness of the proposed hierarchical matrices, we conduct some numerical experiments. The numerical results show that they outperform their uniform admissibility counterparts in terms of storage and matrix-vector product cost, demonstrating the effectiveness of the proposed hybrid formats. We also show that the proposed formats are competitive with the standard (or strong) admissibility-based $\mathcal{H}^2$-matrices. Furthermore, we showcase some interesting applications of the proposed hierarchical matrices.
We present an adaptive bitwidth allocation method for quantizing low-rank factors in hierarchical matrices in order to decrease storage while meeting a given constraint on their reconstruction error. For the bitwidth allocation we formulate a discrete storage minimization problem and show that an optimal assignment can give larger singular values at least as many bits as smaller ones. We further develop a greedy algorithm that meets the specified max-norm quantization error budget. We investigate a hierarchical block partition in which all of the non-diagonal blocks are represented by quantized low-rank factors and only dense blocks remain dense; in this way a global reconstruction bound related to $η_k$ and $R_k$ is obtained, with the former and latter being bounds on the truncation and reconstruction arithmetic errors, respectively. Extensive simulations, carried out with kernel and BEM matrices under block-relative budgets, indicate that the proposed method achieves a reduction in storage of up to $44.3\times$ when compared to storing the complete matrix densely in FP64, and a reduction of up to $5.70\times$ when compared to the FP64 low-rank factors within the same partition, together with measured reconstruction errors that remained below the maximum norm bounds in all the hierarchical tests.
Multiply-with-carry (MWC) generators are very fast pseudorandom number generators whose quality depends on the choice of a multiplier, and is usually assessed by applying the spectral test to an associated multiplicative congruential generator. The special form of the modulus is shown to impose structural limitations on the lattice structure of every MWC generator, which the generalized MWC generators of Goresky and Klapper mitigate; in the tests reported, a scrambler made of a single xor removes them. Since these limitations make the usual worst-case criterion useless, candidates are ranked by a weighted geometric mean of their figures of merit, which does not depend on their normalization, and multipliers are suggested for generators with 64 to 256 bits of state. Along the way, some errors in the literature are corrected. In the end, MWC192X is proposed: a 64-bit generator with 192 bits of state that passes all known tests, and emits an output in half a nanosecond. Analogous 32-bit generators are proposed for GPUs.
Reduced-precision arithmetic offers significant opportunities to improve performance, memory usage, and energy efficiency in numerical applications, provided that numerical accuracy is preserved. This work investigates automated precision tuning through customized floating-point formats with user-defined exponent and significand sizes, enabling the emulation of emerging low-precision formats and the exploration of non-standard precision configurations within a unified mixed-precision framework. The proposed methodology, implemented in the PROMISE precision autotuning tool, combines numerical validation with a systematic search to generate program variants that satisfy user-defined accuracy requirements. To address the computational cost of this exploration, a containerized benchmarking framework supports parallel execution across multiple algorithms and parameter configurations. The approach is evaluated on a suite of numerical programs, including linear solvers and applications from the Rodinia benchmark. Results show that a substantial proportion of variables can be safely reduced to lower precision while preserving accuracy, indicating that standard double precision is often over-provisioned. These findings highlight the potential of automated precision tuning to derive efficient mixed-precision configurations tailored to application-specific accuracy requirements.
Low-rank methods can reduce the memory and computational requirements of deep neural network training in approaches such as GaLore. Randomized range finders offer an attractive alternative to singular value decompositions, particularly when the required rank is determined adaptively from a prescribed approximation tolerance. We introduce a blocked adaptive randomized range finder based on a Householder QR decomposition and an optimized hybrid CPU--GPU implementation. Numerical experiments show that it preserves orthogonality and approximation accuracy in regimes where block Gram--Schmidt without reorthogonalization becomes unstable. The blocked formulation exposes matrix--matrix operations and enables overlap of CPU panel factorization with GPU updates. On an NVIDIA GH200, the overlapped implementation reduces the runtime for the largest tested matrix from 9.91 s on the CPU to 0.407 s. The method provides a stable and efficient building block for low-rank approximation on heterogeneous systems with applications in computational science and engineering.
In this paper, we introduce a practical GPU-enhanced matrix-free first-order method for solving large-scale conic programming problems, which we refer to as PDCS, standing for the Primal-Dual Conic Programming Solver. Problems that it solves include linear programs, second-order cone programs, convex quadratic programs, and exponential cone programs. The method avoids matrix factorizations and uses sparse matrix-vector multiplication as its core operation, which is memory-efficient and well suited for GPU acceleration. It builds on a variant of the primal-dual hybrid gradient method (PDHG) called reflected restarted Halpern PDHG with averaging, which we prove converges on general conic programs, and further incorporates several practical enhancements. Its GPU implementation, cuPDCS, uses customized computational schemes that assign cones of different types and sizes to different levels of the GPU architecture. Numerical experiments demonstrate that cuPDCS is generally more efficient than state-of-the-art commercial solvers and other first-order methods on large-scale conic program applications, including Fisher market equilibrium problems, Lasso regression, and multi-period portfolio optimization. cuPDCS also exhibits better scalability, efficiency, and robustness compared to other first-order methods on the conic program benchmark dataset CBLIB. These advantages are more pronounced in large-scale settings.
Low-precision GPU kernels are commonly tested against a reference using hand-picked absolute and relative tolerances. Legitimate numerical error varies with the input, format and precision plan, so a fixed tolerance trades false alarms for missed faults. We present numoracle, which derives per-element acceptance regions from precision-annotated reference programs. Its calculus tracks shared rounding errors, provenance-ordered products, block-scaled quantizers and calibrated math-library errors. A rounding-interval check detects faulty stores that symmetric bounds can admit. Soundness is conditional on explicit rounding, hardware, domain and implementation assumptions; a fixed-input mode removes the zero-mean rounding assumption. We evaluate 13 Triton operators from Liger-Kernel, Unsloth and vLLM on two GPU generations. The main campaigns produce no false alarm on 33,960 correct-kernel and behavior-preserving-variant runs per GPU. Under input shift, the GPU campaigns record no false alarm while tuned tolerances flag up to 30% of correct runs; final-version reference checks identify inputs without complete finite-range certificates. On the RTX PRO 6000, numoracle kills 116 of 124 effective mutants, versus 59-83 for four fixed baselines. Among 400 tuned tolerance configurations, requiring at most 1% false alarms in every evaluated Liger input group permits at most 31 of 53 main-mutant kills, versus numoracle's 50. Historical replays expose accuracy faults in eight change pairs without rejecting their matched corrected controls; a ninth change concerns rounding parity.
We present cuPDQPx, a GPU-accelerated, factorization-free first-order solver for large-scale convex quadratic programming (QP). Building on the restarted reflected Halpern primal-dual framework of cuPDLPx, we introduce an auxiliary-variable splitting that yields explicit updates without linearizing the quadratic objective or solving inner quadratic subproblems. A preconditioned proximal point interpretation yields both sublinear and linear fixed-point residual convergence guarantees. For practical efficiency, cuPDQPx integrates diagonal preconditioning, step-size selection, adaptive restart, and adaptive primal and quadratic weight updates. cuPDQPx is implemented in C and CUDA with optimized sparse matrix-vector products, kernel fusion, and CUDA Graphs. Numerical experiments on the Maros-Mészáros collection, a convex subset of QPLIB, and synthetic large-scale QPs demonstrate that cuPDQPx has strong numerical performance at both $10^{-4}$ and $10^{-8}$ tolerances.
Antoine Poirot-Bourdain, Alexis Van Straaten, Emmanuel Bacry, Anne-Sophie Jannot, Jean Feydy
This paper introduces survivalGPU, a high-performance Python and R package designed to accelerate survival analysis computations using graphics processing units (GPUs). While the Cox proportional hazards model remains a cornerstone of survival analysis, it can impose substantial computational burdens when applied to large longitudinal data sets. To address this limitation, survivalGPU implements Cox proportional hazards models and weighted cumulative exposure (WCE) models in Python, using PyTorch and KeOps for efficient back-end computation, with interfaces in both Python and R. Numerical validation shows that survivalGPU reproduces the results of reference implementations within the tolerances expected for the selected floating-point precision, with relative errors below 10 -2 in float32 and below 10 -5 in float64. Benchmark experiments demonstrate substantial speedups. For Cox models, coxphGPU() was 5.6x faster than coxph() from the survival package on 1,000,000 patients, and 14.0x faster based on runtimes extrapolated to 1,000 bootstrap replicates on 100,000 patients. Gains were larger for WCE models: wceGPU() was 920.7x faster than WCE() from the WCE package on 100,000 patients, and 1,321.4x faster based on runtimes extrapolated to 1,000 bootstrap replicates. These results show that survivalGPU makes flexible time-dependent survival modeling and bootstrap-based estimation of confidence intervals feasible for large-scale pharmacoepidemiological applications.
Learning physical models from spatiotemporal data requires simulating systems of partial differential equations (PDEs) and differentiating through their solution. Upon spatial discretization, these models give rise to stiff dynamical systems that require repeatedly solving large, sparse linear systems. Traditional scientific computing frameworks are largely CPU-based and require problem-specific tuning to accelerate, scale, and differentiate. In this work, we introduce diff-soup, a differentiable, sparse-oriented, unified PDE suite for GPU-native simulation and inference of large, stiff dynamical systems. Built in JAX, diff-soup provides 1) a fully differentiable design for problems with sparse structure, 2) scaling to large problems through the NVIDIA cuDSS sparse direct solver with model and data multi-GPU parallelism, and 3) end-to-end PDE-constrained optimization, from classical statistical inference to embedded deep learning. On stiff PDE benchmarks, diff-soup outperforms established CPU- and GPU-based frameworks by one to two orders of magnitude in both forward solves and gradients. We demonstrate both a classical Bayesian inference workflow, identifying constitutive relationships, and a modern deep learning workflow, training a neural field for the catalyst activity of a three-dimensional packed bed reactor. We provide diff-soup as an open-source repository to facilitate research in computational science and engineering.
The Low Order-Value Optimization (LOVO) problem involves minimizing the minimum among a finite number of function values within a feasible set. LOVO has several practical applications such as robust parameter estimation, protein alignment, portfolio optimization, among others. In this work, we are interested in the constrained nonlinear optimization LOVO problem of minimizing the minimum between a finite number of function values subject to a nonempty closed convex set where each function is a black-box and continuously differentiable, but the derivatives are not available. We develop the first derivative-free trust-region algorithm for constrained LOVO problems with convergence to weakly critical points. Under suitable conditions, we establish the global convergence of the algorithm and also its worst-case iteration complexity analysis. An initial open-source implementation using only linear interpolation models is developed. Extensive numerical experiments and comparison with existing alternatives show the properties and the efficiency of the proposed approach when solving LOVO problems.
Automatic differentiation (the automatic computation of numerical function derivatives) has many applications, especially in scientific computing and machine learning. Several approaches to automate the calculation of derivatives have been presented. Common challenges facing those approaches include (1) having to define derivative rules for different classes of functions, (2) having to manipulate program symbols to understand the semantics of the function being differentiated, and construct the derivative function, and (3), more importantly, ensuring that the calculated derivative is correct. In this paper, we present an approach for automatic \emph{certified} differentiation using the Lean proof assistant. By certified, we mean a correctness certificate is computed together with the derivative function and a domain on which the derivative is guaranteed to be valid. Our approach leverages existing derivative theorems included in the Mathlib Lean library instead of having to define our own rules. We also utilize Lean metaprogramming to manipulate the symbolic structure of the input function, and to construct the derivative function. The correctness certificate is computed step-wise as we apply Mathlib theorems. To make the search efficient, we index the derivative theorems by the function they differentiate. We evaluate our approach on 500 randomly generated functions: with the index, all of them are differentiated, with a median time of 0.18\,s per function, 8.6 times faster than without the index.
Mathematical software has particular needs that are not all covered by the existing literature on software engineering or algorithm design. This paper looks at some of these aspects, in the concrete context of the computer algebra system GAP. It also studies the potential impact AI systems could have on the development process for mathematical software.
We describe a new statistical test for pseudorandom number generators (PRNGs). Our test can find bias induced by dependencies among the Hamming weights of the outputs of a PRNG, even for PRNGs that pass state-of-the-art tests of the same kind from the literature, and in particular for generators based on $\mathbf F_2$-linear transformations such as the dSFMT, xoroshiro1024+, and WELL512a.
Recovering the governing equations of a dynamical system from data is a central problem in data-driven science, yet the available Python tooling forces a compromise among automation, statistical rigor, and reliability. We present py-ARGOS, a Python package implementing two automatic algorithms, bootstrap ARGOS and Bayesian-ARGOS, behind a uniform interface. Both are statistically principled throughout: automated model selection reduces the candidate library, each coefficient is reported with a bootstrap confidence or Bayesian credible interval that quantifies its uncertainty, and a term is retained when this interval excludes zero. Built-in diagnostics additionally flag unreliable fits. On chaotic benchmark systems, the two algorithms improve on the SINDy family's identification performance.