pith. sign in

cs.MS

Mathematical Software

Roughly includes material in ACM Subject Class G.4.

0
math.NA 2026-07-01

Any-degree Qr factorization triangularizes polynomial systems

by Dileepan Joseph

Ricci-Notation Tensor Framework for Numerical Algebraic Geometry via Any-Degree Unitary-Triangular Factorization

The unitary factor preserves the original solution set for zero- and positive-dimensional cases, extending numerical methods beyond linear a

Figure from the paper full image
abstract click to expand
The unitary-triangular (QR) factorization of linear algebra may be used to robustly and efficiently solve a linear system. Toward a comparable numerical method to solve a polynomial system of higher degree, this paper proposes an any-degree unitary-triangular (Qr) factorization, which for a degree-one system reduces to the QR factorization. The work develops a tensor framework, i.e., codesigned algebra and software, where polynomial system coefficients are represented by a vector-shaped sparse tensor, a multidimensional array whose number of Ricci-notation indices, called the tensor degree, equals the highest monomial degree of the system. With the proposed Qr factorization, the coefficient tensor decomposes into a product of unitary and triangular factors that, in general, also have Ricci-notation indices and sparse entries. The unitary factor defines a unitary transform, a generalization of the related linear algebra concept to tensor algebra, that can triangularize a polynomial system while preserving its solution set, whether zero- or positive-dimensional. The work extends the author's Ricci-notation tensor framework, providing new algebra and new software to model, construct, and factorize polynomial systems in this manner. After applying the approach to numerically triangularize two zero-dimensional systems, chosen for educational value, results are compared to the Gr\"obner-basis (GB) method for triangularizing polynomial systems symbolically. One problem is of degree three, with three equations and unknowns, and the other of degree two, with four equations and unknowns. Although it resembles GB triangularization, the proposed Qr factorization has a completely different pedigree associated with numerical methods.
0
0
cs.MS 2026-06-30

GPU kernel cuts Trotter time 7.6x for 2D Ising tensor networks

by Paolo D'Alberto

GPU-First Heisenberg-Picture Tensor Network Dynamics for the 2D Transverse-Field Ising Model

Zero-malloc workspace and device-only permutations remove host bottlenecks during Heisenberg-picture contraction.

Figure from the paper full image
abstract click to expand
We present CppSim, a C++/GPU 2D Ising simulator for Heisenberg-picture tensor network time evolution on GPUs. The key computational contributions are: first, a zero-malloc GPU workspace that pre-allocates all buffers at startup; second, a custom GPU tensor permutation kernel replacing host-side index shuffling with a pure device-to-device operation, yielding a 7.6x trotter speedup; third, a hybrid QR strategy selecting Cholesky-QR for tall-skinny matrices and Householder-QR otherwise; fourth, adaptive Belief Propagation with log-space Bethe partition function evaluation and explicit sign tracking.
0
0
math.OC 2026-06-30

Library solves SPPRC 1.3-2.35 times faster than PathWyse

by Simon Spoorendonk

texttt{bucket-graph-spprc}: an extensible C++ library for the shortest path problem with resource constraints

bgspprc uses compile-time resources and bucket-graph labelling for faster pricing in vehicle routing solvers.

Figure from the paper full image
abstract click to expand
We present $\texttt{bucket-graph-spprc}$ ($\texttt{bgspprc}$ for short), an open-source, header-only C++23 library for the shortest path problem with resource constraints (SPPRC), the pricing subproblem at the heart of branch-cut-and-price for vehicle routing and related problems. The library implements the bucket-graph labelling algorithm of Sadykov, Uchoa and Pessoa (2021), with bidirectional labelling, across-arc concatenation, bucket fixing and arc elimination, and a structure-of-arrays label store with SIMD-accelerated dominance. Its central design feature is a compile-time resource concept: a new SPPRC variant is added by implementing a fixed seven-function interface, and resources compose into a label state with no runtime dispatch, the state layout fixed at compile time. Five resources ship built in: time/capacity, ng-path elementarity relaxation, rank-1 cuts, cumulative cost, and pickup-and-delivery. In a reproducible, head-to-head comparison on shared public instances at an identical bound, $\texttt{bgspprc}$ outperforms PathWyse (Salani, Basso and Giuffrida, 2024), the main open-source comparator, by $1.3\times$--$2.35\times$ in shifted geometric mean (and by $1.3\times$--$2.3\times$ even when itself run single-threaded), and runs within $1.9\times$--$2.4\times$ of parallel pull labelling (Petersen and Spoorendonk, 2025), a different labelling technique for the same problem. The library, benchmark scripts, and pinned instances are publicly available.
0
0
cs.MS 2026-06-29

Revised scaling fixes scale invariance in Ozaki fast mode

by Shota Kawakami, Daisuke Takahashi

Improved Scaling for Fast Mode of Ozaki Scheme II

New formula from CRT uniqueness condition guarantees validity for any input scaling at no extra cost and matches accurate-mode accuracy at f

Figure from the paper full image
abstract click to expand
Ozaki scheme II emulates high-precision matrix multiplication using low-precision integer matrix operations based on the Chinese remainder theorem (CRT). It first scales the high-precision matrices to convert them into integer matrices. For this scaling step, Ozaki scheme II provides two modes: accurate mode, which uses INT8 matrix multiplication to estimate scaling factors, and fast mode, which applies the Cauchy--Schwarz inequality at lower computational cost. We show that the existing formula lacks scale invariance; multiplying the input matrices by a constant changes the effective bit width of the integer matrices in the scaling step, causing accuracy degradation or CRT recovery failure. To address this, we propose a revised scaling formula derived from the CRT uniqueness condition via the Cauchy--Schwarz inequality. The proposed formula is scale-invariant by construction, guarantees that the CRT uniqueness condition is always satisfied, and introduces no additional overhead over the original fast mode. Experiments on an NVIDIA GH200 GPU show that the proposed method achieves accuracy comparable to that of accurate mode while maintaining throughput comparable to that of fast mode. In the accuracy--throughput trade-off, the proposed method overcomes the accuracy limitation of fast mode and the throughput constraint of accurate mode, offering a superior accuracy and performance.
0
0
cs.MS 2026-06-29

Package computes alpha-stable densities more reliably

by Wael Tabbara, Sharafeddine Sharafeddine +3 more

Density Functions and Random Number Generators of α-Stable Distributions

AUB-HTP adds multivariate sampling and widens the usable parameter range for these heavy-tailed distributions.

Figure from the paper full image
abstract click to expand
Heavy-tailed distributions are increasingly found to better fit empirical data in engineering, finance, physics, network science, and related fields. Among them, $\alpha$-stable distributions play a central role being limiting laws in the generalized central limit theorem: they are expected to be exceptionally good models whenever sums of multiple independent heavy-tailed sources are at play. Despite their theoretical importance, their practical use remains challenging: $\alpha$-stable probability densities generally do not have closed-form expressions, and numerical evaluation and random variate generation can be difficult, especially in the multivariate setting. This paper presents AUB-HTP, a Python package for numerical computation and simulation of $\alpha$-stable distributions. The package provides scalar density evaluation using several complementary methods, including Zolotarev-type integral representations, series formulas, and numerical inversion of characteristic functions. It also provides random variate generation for scalar and multivariate $\alpha$-stable distributions, with support for flexible spectral measures through LePage series representations. Numerical experiments demonstrate that AUB-HTP improves the accuracy, stability, and parameter coverage of existing tools for scalar density computation, while adding new capabilities for multivariate simulation. The package is designed to support reproducible computational work involving heavy-tailed models across a broad range of scientific applications.
0
0
physics.comp-ph 2026-06-26

Python library designs constrained discrete filters automatically

by Z. Nikolaou, P. Domingo +2 more

pyDOF: a Python library for the design of discrete forward and inverse filters

Users set monotonicity, positivity and other rules; coefficients are written to plain text usable in any code.

Figure from the paper full image
abstract click to expand
In this work, we present pyDOF, a Python-based software library which provides a domain-specific framework for the design of symmetric, physical-space, forward as well as inverse discrete filters. pyDOF is based on a constrained optimisation framework developed in our previous work [1, 2]. This framework allows the user to impose a wide range of constraints on the discrete filter transfer-function such as monotonicity, positivity, value-fixing, gradient-smoothing etc. amongst many others. pyDOF additionally includes an adaptive filter stencil selection option, and a van Cittert-based inverse-filter design with a user-controlled reconstruction order. The filter coefficients are computed automatically, and saved to a plain text file which can be readily parsed by any programming language. pyDOF can be used to design a wide range of low-pass, high-pass, multi band-pass/band-stop etc. discrete filters. In addition, due to its generality and abstraction, pyDOF can be used to design specific filters for user-defined target filter transfer functions. Although developed primarily for application to computational fluid dynamics simulations, pyDOF can be used to design discrete filters for a wide range of signal processing applications.
0
0
cs.DC 2026-06-25

Fused kernels reach 83% of INT8 peak for emulated high-precision GEMM

by Denghui Lu, Alexander Maeder +2 more

EmuGEMM: Fused Tensor Core Kernels for Precision Emulation in Matrix Multiplication

By keeping Ozaki intermediates on-chip, EmuGEMM beats cuBLAS TF32 by up to 1.7x on Hopper and Blackwell at matching accuracy.

Figure from the paper full image
abstract click to expand
Modern GPUs devote an increasing silicon budget to low-precision matrix-multiplication units, widening the precision-throughput gap for scientific computing workloads. Ozaki Schemes I and II offer an alternative by reconstructing high-precision general matrix multiplication (GEMM) from low-precision operations, yet existing implementations leave substantial performance untapped. In particular, intermediate results are repeatedly materialized in global memory, making data movement the dominant bottleneck. We present EmuGEMM, fused integer Tensor Core kernels for NVIDIA Hopper and Blackwell GPUs that eliminate redundant memory round-trips in both Ozaki schemes. Using Scheme I, EmuGEMM sustains up to 1,639 Top/s on Hopper (83% of INT8 peak) and 3,654 Top/s on Blackwell (81%). For large matrices, EmuGEMM surpasses cuBLAS TF32 throughput by up to 1.4x on Hopper and 1.7x on Blackwell, at comparable accuracy. Using Scheme II, EmuGEMM extends to complex arithmetic and outperforms cuBLAS ZGEMM by up to 2.3x on Hopper and 5.5x on Blackwell.
0
0
math.OC 2026-06-25

Zoutendijk theorem formalized in Lean for line search

by Yiyang Zhang, Kenneth W. Shum

Formalization of Line Search Methods by Lean

Machine-checked proofs cover Armijo, Wolfe conditions and global convergence for gradient-based methods

abstract click to expand
This paper presents a formalization of line search methods in the Lean 4 theorem prover. Our goal is to advance machine verification of nonlinear optimization theory by translating standard textbook definitions and convergence arguments into rigorous Lean code. We formalize fundamental notions related to gradient descent and descent directions, adaptive step-size selection via backtracking line search, and several classical line search criteria, including the Armijo, Goldstein, and Wolfe conditions, as well as nonmonotone variants. We further formalize a key convergence result, namely the Zoutendijk theorem, which plays a central role in the global convergence analysis of gradient-based iterative methods. By providing machine-checkable definitions and proofs for line search theory, this work complements existing formalizations of first-order optimization methods and establishes a foundation for the verified development of more advanced algorithms in nonlinear programming.
0
0
cs.MS 2026-06-24

Carleman truncations converge geometrically over certified time horizons

by Michele Boreale, Luisa Collodi

New convergence results for Carleman linearization

Explicit degree-aware bounds from Dyson-Duhamel expansion give finite-degree estimates and improve on prior results for test ODEs.

Figure from the paper full image
abstract click to expand
We prove new error bounds for finite Carleman truncations of polynomial ordinary differential equations. The analysis works directly in the original monomial basis and for selected observables, such as state coordinates. Using a Dyson--Duhamel expansion, we separate the degree-preserving linear part from the degree-raising nonlinear part and track how truncation errors can propagate back to the observable. The resulting bounds are degree-aware and retain logarithmic-norm information from the original linear dynamics. We obtain explicit finite-degree estimates and geometric convergence over certified time horizons. Comparisons with existing bounds, in particular those of Forets--Pouly, are given on the Stuart--Landau and Van der Pol systems.
0
0
cs.DC 2026-06-23

Struct splitting trims SPH GPU packing time by 20-40%

by Mladen Ivkovic, Abouzied M.A.Nasar +5 more

Memory Layouts for GPU-Data Transfer Buffering in SPH

Access-pattern decomposition of particle data lowers total offloading overhead by 12-25% as transfers dominate runtime.

Figure from the paper full image
abstract click to expand
The rise in GPU compute speed has outpaced improvements in host-to-device memory transfer speeds, despite the advent of shared-memory superchips. Consequently, memory transfer times now constitute an increasingly large fraction of total time-to-solution, compelling developers to compress GPU kernel input and output data into compact, minimal formats prior to GPU-offloading. This complements existing work on GPU- and compute-friendly data arrangements. We study a Smoothed Particle Hydrodynamics solver and propose memory layout strategies for host-side particle data that are particularly well-suited to GPU-offloading. Specifically, we advocate splitting classic array-of-struct data structures into a split array-of-struct arrangement, in which each logical struct decomposes into substructs determined by kernel read/write access patterns and attribute types. Splitting a monolithic particle struct into several bespoke, finer-grained structs can reduce the time required to pack data to and from buffers by ~20% - 40%, lowering total time spent on GPU-offloading by ~12% - 25%.
0
0
cs.CG 2026-06-22

Iterative method produces arc-length splines that hit all points

by Dafna K. Matsegora, Stephen M. Watt

Arc-Length Parameterized Interpolating Splines

Previous techniques either miss the points or use non-arc-length parameterization; this keeps both.

Figure from the paper full image
abstract click to expand
We present an iterative algorithm to compute an arc-length parameterized spline interpolating a set of points. This differs from other methods where the computed spline either does not interpolate the original points or the parameterization is not the arc-length of the returned curves. Our method is applicable in any dimension $D \ge 2$, and we illustrate it with numerical results for plane curves.
0
0
cs.MS 2026-06-22

Asynchronous Taylor solver scales linearly for large DDEs

by Avinash Malik

An Asynchronous multi-rate Taylor method for Delay Differential Equations

Local clocks and static buffers advance only active sub-graphs while eliminating dynamic memory allocation and interpolation.

abstract click to expand
The numerical simulation of high-dimensional, multi-rate Delay Differential Equations (DDEs) is fundamentally bottlenecked by synchronous time-stepping and the dynamic memory allocation required for continuous history tracking. In this paper, we introduce the Asynchronous Adaptive Taylor Solver (AATS), an event-driven integration framework designed to overcome these high-performance computing limitations. By assigning independent local clocks to individual coordinates and advancing them using high-order Taylor polynomials generated via compile-time Automatic Differentiation, AATS restricts computational work to actively evolving sub-graphs. To eliminate the severe memory overhead endemic to traditional DDE solvers, AATS utilizes statically allocated circular buffers to store polynomial segments, achieving interpolation-free continuous dense-output evaluation with a verified zero-allocation runtime memory footprint. Alongside this software architecture, we establish a novel continuous proof of convergence for asynchronous Taylor expansions and formally prove that the framework's algorithmic complexity scales linearly (O(N)). Extensive benchmarks against state-of-the-art synchronous solvers (Julia SciML) validate these theoretical bounds. On large-scale benchmarks (upto $N = 10000$ coordinates) AATS fundamentally minimizes the constant factor of algorithmic work by avoiding redundant evaluations, delivering empirically consistent with O(N) execution scaling and significant wall-clock speedups.
0
0
cs.MS 2026-06-19

MLIR dialect compiles finite element ops for CPU and GPU

by Edward Erasmie-Jones, Giacomo Castiglioni +1 more

NektarIR: A Domain-Specific Compiler for High-Order Finite Element Operations on Heterogeneous Hardware

High-level abstractions for spectral/hp operators lower progressively with domain optimizations at each stage.

Figure from the paper full image
abstract click to expand
Modern high performance computing (HPC) applications must target heterogeneous hardware. This requires significant work to ensure domain specific implementations translate to highly performant kernels across a range hardware types and vendors, each requiring bespoke optimization to make use of the specific target architecture. Through the development of a domain specific compiler built with the multi-level intermediate representations (MLIR) project, one can express a high-level, close to the specific domain, abstraction that is progressively lowered to a low, close to metal, abstraction. At each intermediate representation (IR), appropriate optimizations can be applied without costly analysis due to the knowledge embedded in the domain specific IRs. We apply this method to the construction of discrete differential operators for use in spectral/hp element method solvers for computational fluid dynamics (CFD). Here, the performance is driven by a small set of common finite element operators that are composed to create kernels for the discrete differential operators used to solve weak partial differential equations. We create our own MLIR dialect to represent these operators and implement a bespoke lowering pipeline to facilitate the just-in-time compilation of these kernels for both CPU and GPU architecture and illustrate performance comparisons with the Nektar++ spectral/hp element framework.
0
0
cs.MS 2026-06-19

Multiword sums of doubles pair with parallel hardware

by Jan Verschelde

Multiword Arithmetic and Parallel Computing

Register storage of unevaluated 64-bit sums supports GPU tasking for higher precision without custom libraries.

Figure from the paper full image
abstract click to expand
In many applications, the precision by the available hardware arithmetic is insufficient to guarantee accurate results. Multiword arithmetic is a special type of multiprecision arithmetic where a multiple double is an unevaluated sum of 64-bit doubles, or where a multiple integer is an unevaluated sum of 64-bit integers. Parallel computing is applied to compensate for the cost overhead of multiword arithmetic. This type of arithmetic exploits naturally the optimized hardware, allows for efficient type conversions, memory layouts, all favorable for parallel computing. For example, storing a multiword in registers rather than arrays is beneficial to parallel computing by tasking and acceleration by graphics processing units. Code for multiword arithmetic is available in the software PHCpack, written mainly in Ada, publicly available at github, and as an Alire crate, released under the GNU GPL v3.0 license.
0
0
cs.MS 2026-06-19

Compiler permits any parallel mix of sparse tensor operations

by Adrian Gushin, Sang Yoon Kim +1 more

WingSpan: Concurrency and Dependence for Sparse and Structured Tensor Compilers

WingSpan reaches hand-tuned SpGEMM speed and uses dependence analysis to keep compositions safe.

Figure from the paper full image
abstract click to expand
Sparse tensors represent data that is mostly zero or some other compressible fill pattern. Such datasets can be massive, so optimized tensor algebra libraries and compilers have been developed to exploit these patterns to improve performance. Existing systems, however, frequently lack support for parallelism, especially when outputs are sparse or multiple inputs are sparse. We propose WingSpan, a sparse tensor language enabling unrestricted parallel programming. WingSpan supports arbitrary composition of parallel loops and data structures, matching or exceeding the performance of hand-optimized parallel routines on critical kernels such as SpGEMM. We also introduce a dependence theory for the safety of parallel programs involving sparse tensors and structures beyond sparsity.
0
0
math.NA 2026-06-19

Structured coarse space removes communication from Poisson solver

by Thilina Ratnayaka, Paul Fischer +1 more

Coarse Solvers for Exascale Solution of Poisson Problems

Non-nested global space enables communication-free interpolation inside two-level Schwarz method for pMG pressure Poisson solves

Figure from the paper full image
abstract click to expand
We present a two-level Schwarz method as an alternative to Algebraic Multigrid method(AMG) used as the last level (coarse) solver of the p-multigrid pMG preconditioner for pressure Poisson equation resulting from Spectral/Finite element descretization of incompressible Navier-Stokes equation. Proposed Schwarz method consits of a local problem in the original pMG coarse space and a global coarse problem. Main contribution of the paper is a novel, structured and a non-nested coarse space for the global coarse problem. Structured nature of the proposed global coarse space enable communication-free interpolation between the original p-multgrid coarse space and the global coarse problem. We demonstrate the effectiveness of the proposed method compared to the state of the art AMG solver BoomerAMG by a series of experiments performed using Nek5000/RS, a suite of highly scalable incompressible Navier-Stokes solvers, on Summit/Frontier supercomputers at Oak Ridge Leadership Computing Facility.
0
0
cs.MS 2026-06-19

Unified interface reduces solver binding effort for optimization tasks

by Dmitry I. Kabanov, Stephan Rave +1 more

Software package MaRDI Open Interfaces for improved interoperability in numerical optimization

MaRDI Open Interfaces standardizes connections to different solvers, freeing researchers to focus on experiments like PINN training for PDEs

Figure from the paper full image
abstract click to expand
To address the challenges of interoperability in computational science, we present the latest updates to the software package MaRDI Open Interfaces. This software package aims to decrease the time and coding/testing efforts spent by computational scientists on tasks such as writing bindings to numerical solvers and adapting experiment codes to the varying interfaces of solvers for the same problem type (e.g., for benchmarking, which solver is better). By streamlining these tasks, this software package helps researchers focus on the actual essence of their computational projects. Here, we demonstrate a recently developed interface for nonlinear optimization and illustrate how it can be applied for computational experiments with optimization problems. As an example of such problem, we consider training of physics-informed neural networks to predict the solutions of viscous Burgers' equation.
0
0
cs.MS 2026-06-19

Unified interfaces speed solver switches in nonlinear optimization

by Dmitry I. Kabanov, Stephan Rave +1 more

MaRDI Open Interfaces for Interoperable Nonlinear Optimization

Automated marshalling cuts code changes when testing multiple solvers across languages.

Figure from the paper full image
abstract click to expand
MaRDI Open Interfaces is a software package that aims to improve interoperability in scientific computing, particularly, for nonlinear optimization. To this end, this package holds two main characteristics. First, it provides unified interfaces for typical numerical problems to help switching between solvers for the same problem type. Second, it automates data marshalling between programming languages. Hence, computational scientists can conduct experiments faster by using the package, with fewer code-modification and testing efforts. In this work we describe the general structure of the software package and show examples with the interface for nonlinear optimization.
0
0
math.CV 2026-06-19

Lean 4 verifies uniqueness of Möbius maps on three points

by Fubin Yan, Kenneth W. Shum

Formalizing Extended Complex Numbers, Mobius Transformations, and Cross Ratio in Lean 4

Machine-checked proofs confirm that any three distinct points determine a unique Möbius transformation and that the cross ratio stays fixed.

abstract click to expand
The extended complex plane is a fundamental object in complex analysis, hyperbolic geometry, and mathematical physics. Its geometry is governed by M\"obius transformations, with the cross ratio serving as a central invariant. We present a formalization of these concepts in the Lean4 theorem prover. The extended complex plane is represented using Mathlib's Option type over $\mathbb{C}$, where the additional element represents the point at infinity. On this foundation, we define M\"obius transformations, their action on the extended complex plane, and the cross ratio. We formalize several basic properties of M\"obius transformations, including their group structure, and identify them with a projective general linear group. We also prove the uniqueness of a M\"obius transformation mapping any three distinct points to any other three distinct points, and the invariance of the cross ratio. All proofs are machine-checked in Lean 4. The complete development comprises approximately 6,000 lines of Lean code, including about 40 definitions and 150 lemmas and theorems. This work provides a verified foundation for future formalizations of conformal geometry, hyperbolic models, modular forms, and applications in mathematical physics.
0
0
cs.MS 2026-06-18

Rust matches Eigen and PSBLAS on sparse matrix kernels

by Luca Lombardo, Fabio Durastante

Evaluating Rust for Sparse Matrix Kernels in Scientific Computing

Benchmarks across SpMV, Lanczos, and exponential workloads show native Rust code tracks state-of-the-art CSC performance.

abstract click to expand
Sparse matrix kernels form the computational backbone of scientific computing, traditionally relying on C/C++ and Fortran implementations that prioritize performance over memory safety. This work evaluates Rust as a systems-level alternative for sparse linear algebra by implementing and benchmarking three core workloads: sparse matrix-vector multiplication (SpMV), Lanczos-based Krylov methods, and matrix-exponential evaluation. We compare native Rust code against established baselines (Intel oneMKL, Eigen, PETSc, and PSBLAS) across a suite of representative matrices. Our results show that Rust's sparse kernels achieve performance comparable to Eigen and PSBLAS, tracking the state-of-the-art for CSC formats, while trailing PETSc's advanced blocked CSR optimizations. By analyzing compile-time monomorphization, SIMD vectorization, and FFI boundaries, we assess the practical impact of Rust's safety model and ecosystem readiness. The study provides concrete, evidence-based guidance for modernizing high-performance numerical software stacks.
0
0
stat.ME 2026-06-18

Neural matching estimates individual treatment effects in time series

by Maha Ouali, Badih Ghattas +3 more

Balanced Twins: Causal Inference on Time Series with Hidden Confounding

By learning latent representations and propensity scores together, the method matches units to recover counterfactuals under hidden bias and

Figure from the paper full image
abstract click to expand
Accurately estimating treatment effects in time series is essential for evaluating interventions in real-world applications, especially when treatment assignment is biased by unobserved factors. In many practical settings, interventions are adopted at different times across individuals, leading to staggered treatment exposure and heterogeneous pre-treatment histories. In such cases, aggregating outcome trajectories across treated units is ill-defined, making individual treatment effect (ITE) estimation a prerequisite for reliable causal inference. We therefore study the problem of estimating the average treatment effect for the treated (ATT) by first recovering individual-level counterfactuals. We introduce a neural framework that learns simultaneously low-dimensional latent representations of individual time series and propensity scores. These estimates are then used to approximate the individual treatment effects through a flexible matching procedure that avoids classical convexity constraints commonly used in synthetic control methods. By operating at the individual level, our approach naturally accommodates staggered interventions and improves counterfactual estimation under latent bias, without relying on explicit temporal modeling assumptions. We illustrate our approach on both real-world energy consumption data and clinical time series, including high-frequency electricity demand-response programs and semi-synthetic data for individuals in intensive care unit (ICU), where hidden confounding, staggered treatment adoption, and non-stationary dynamics are prevalent.
0
0
cs.MS 2026-06-17

L2 error stays good while PDE simulator diagnostics worsen

by Lennon J. Shikhman

A Diagnostic Software Suite for Auditing Learned PDE Simulators

Suite of structural checks on five benchmarks shows single error score misses coherence failures

abstract click to expand
Learned PDE simulators are increasingly used as low-cost replacements for expensive numerical solvers, but standard relative $L^2$ error does not determine whether a learned model behaves as a coherent numerical time propagator. This paper presents a diagnostic software suite for auditing learned PDE simulators as approximate evolution operators. The suite provides architecture-independent, post hoc diagnostics for relative state error, semigroup consistency, finite-difference generator discrepancy, energy behavior, integral balance, admissibility constraints, perturbation response, and scaling-law consistency. The software is designed around a minimal contract: reference trajectories, a learned propagator or saved predictions, equation metadata, and a diagnostic configuration specifying which structures are meaningful for the problem under study. We validate the suite on five benchmark PDE tasks: two-dimensional incompressible Navier-Stokes, shallow-water dynamics, active matter, three-dimensional compressible Navier-Stokes, and three-dimensional magnetohydrodynamics, using FNO, DeepONet, U-Net, and ResNet-style surrogate models together with controlled underfit and oversmoothed variants. The validation study shows that relative $L^2$ error can remain moderate, or even improve, while structural diagnostics deteriorate substantially. The package therefore supports software-level auditing of learned PDE simulators by reporting an interpretable diagnostic panel rather than collapsing model behavior into a single state-error score.
0
0
cs.PL 2026-06-15

Python package runs models directly on Apple Neural Engine

by Spencer H. Bryngelson

ANEForge: Python for direct computation on the Apple Neural Engine

Compiles lazy tensor graphs to native ANE programs that finish ResNet-18 forward passes in 0.33 ms and support training.

Figure from the paper full image
abstract click to expand
ANEForge is a Python package that programs the Apple Neural Engine (ANE), the fixed-function neural accelerator on every recent Apple device, directly and without CoreML. In production the engine is reachable only through CoreML, which treats it as a scheduling option: no configuration requires the ANE, and a model can silently run on the CPU or GPU instead. ANEForge compiles a lazy tensor graph, built from 58 fused operators and 19 native bridge operators, into a single ANE program. The program is dispatched through the same ANE daemon and kernel-driver stack as Apple's internal framework. Beyond inference, the package reaches the engine's native fused attention, streams int8, int4, and sparse weights, keeps decoder and optimizer state resident across steps, and runs the forward pass, backward pass, and optimizer update of training on the engine. A small fused program completes a call in about 90us, near the engine's 70us per-program dispatch floor, and a pretrained ResNet-18 forward runs end-to-end in 0.33ms. ResNet-18, a sentence encoder, and a Vision Transformer run end-to-end against framework references, and a Stable Diffusion U-Net validates its forward pass. ANEForge targets Apple Silicon under macOS 14 and later. Each release is verified against a recorded macOS and ANE-compiler version.
0
0
cs.AR 2026-06-09

83-format catalog supplies bit-exact packs for FP8 and BF16

by Dmitrii Vasilev

An 83-Format Numeric Catalog with Bit-Exact Conformance Vectors: A Vendor-Neutral Reference for FP8, BF16, MXFP4, and Microscaling Formats

JSON documents and IEEE P3109 cross-walk give engineers a shared reference to diagnose numeric divergences across accelerators.

abstract click to expand
Numeric format proliferation in machine learning hardware -- FP8 (E4M3 and E5M2), BF16, MXFP4, microscaling block formats, and dozens of research variants -- has outpaced the availability of vendor-neutral, bit-exact reference material. Engineers porting models across accelerators encounter silent divergences that are difficult to diagnose without a shared ruler. This paper describes a catalog of 83 numeric formats spanning 13 families, a suite of six bit-exact conformance packs covering GF16, MXFP4 element, BF16, FP8 E4M3, FP8 E5M2, and E8M0 block scale, and an IEEE P3109 v3.2.0 cross-walk that maps each pack to its corresponding standards-track configured format. Each pack is a self-contained JSON document with a SHA-256 fingerprint, a shared row schema, and an anchor vector that encodes 3.0 -- the identity phi^2 + 1/phi^2 = 3 -- as a cross-pack sanity check. Packs are cross-validated against ml_dtypes 0.5.4 (Google/JAX); any divergence is documented explicitly and interpreted as a spec-permitted interpretation gap rather than hidden. The work is framed as registry filling: it does not propose new formats, make model-accuracy claims, or assert superiority over any vendor's implementation. All artifacts are publicly available at https://github.com/gHashTag/t27 under an open license.
1 0
0
cs.MS 2026-06-09

Library wraps AmgX for differentiable GPU AMG solves in JAX

by Yi Liu, Xiantao Fan +1 more

JAX-AMG: A GPU-Accelerated Differentiable Sparse Linear Solver Library for JAX

The wrapper supplies algebraic multigrid and Krylov methods as native JAX primitives that support JIT, adjoint AD, batching and MPI with set

Figure from the paper full image
abstract click to expand
Sparse linear systems from PDE discretizations are central to scientific computing, yet no existing JAX-ecosystem solver simultaneously provides GPU-accelerated algebraic multigrid (AMG), automatic differentiation (AD), and distributed multi-GPU execution. JAX-AMG fills this gap by wrapping the Nvidia AmgX solver suite as a native JAX primitive, exposing AMG and Krylov methods with configurable preconditioners through a unified interface compatible with JIT compilation, reverse-mode AD via adjoint methods, batched solves, and MPI-based distributed execution. Solver caching amortizes setup costs across repeated solves, making JAX-AMG practical for PDE-constrained optimization and inverse problems. The result is a robust, scalable sparse linear algebra layer that integrates seamlessly into differentiable simulation and scientific machine learning pipelines.
0
0
q-bio.QM 2026-06-08

Open code reproduces Mader model lactate outputs

by Katharina Dunst, Vincent Scharf +2 more

MetaboliSim: a Python implementation of the Mader model for dynamic and steady-state simulation of muscular energy metabolism

Dynamic and steady-state Python versions match published references and give consistent MLSS estimates.

Figure from the paper full image
abstract click to expand
The Mader model is the most widely used mathematical framework for muscular energy metabolism in German-language sport science, underpinning lactate diagnostics, maximal lactate steady state (MLSS) estimation and training prescription. Despite decades of use, neither its dynamic ODE formulation nor its steady-state equations have been available as open code, leaving results based on the model impossible to reproduce independently. We close this gap with MetaboliSim, an open-source Python implementation of both formulations: a dynamic model that integrates the five-variable ODE system (phosphate potential, $\dot{V}\mathrm{O}_2$, muscle and blood lactate, and glycogen) with a fourth-order Runge-Kutta scheme, and a steady-state model that computes MLSS power and the lactate-power relationship in one- and two-compartment variants. We verified implementation correctness against published reference values and assessed physiological plausibility across constant-load, step-test, sprint and running protocols. The implementation reproduces the published reference output within stated tolerances and remains numerically stable throughout (halving the time step changes blood lactate by less than 0.01 mmol/L), with both formulations yielding congruent MLSS estimates. Key physiological behaviour ($\dot{V}\mathrm{O}_2$ on-kinetics, lactate accumulation, PCr dynamics and the sub/supra-MLSS separation) emerges directly from the model equations without protocol-specific tuning, and a sensitivity analysis shows MLSS power varying approximately linearly with $\dot{V}\mathrm{O}_{2\max}$ and nonlinearly with $\dot{V}\mathrm{La}_{\max}$. As the first openly available implementation of the complete Mader model (AGPL-3.0), MetaboliSim lets independent groups reproduce, verify and build on published model-based results. Source code: https://codeberg.org/3phos/metabolisim; Platform: https://metabolisim.org
0
0
cs.MS 2026-06-08

Autotuning finds lower-precision configs for most variables

by Xinye Chen, Thibault Hilaire +1 more

Floating-point autotuning with customized precisions

A systematic search plus validation produces mixed-precision variants that meet accuracy targets while cutting memory and energy use in nume

Figure from the paper full image
abstract click to expand
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.
0
0
stat.CO 2026-06-08

R app merges kappa planning and analysis in one interface

by Zheng Gai, Li Xincheng +2 more

CATEKAPPA: An R Shiny Application for Design and Analysis of Consistency Tests Based on the Kappa Statistic for Categorical Responses

CATEKAPPA lets users without programming skills design consistency tests and interpret categorical agreement via a web app.

Figure from the paper full image
abstract click to expand
The kappa statistic is the most widely used measure of inter-rater agreement for categorical data. Despite its popularity, applied researchers often encounter two major hurdles: (i) determining the sample size required to achieve a desired level of agreement with given power, and (ii) computing appropriate kappa coefficients with proper interpretation. Existing R packages such as irr and kappaSize provide these functionalities but require programming skills and lack an integrated, user-friendly interface. We present CATEKAPPA, an R package that bridges this gap by combining sample size planning (via kappaSize) and agreement analysis (via irr) into a single Shiny-based web application. The package supports Cohen's kappa for two raters, Fleiss' kappa for three or more raters, and Light's kappa, and provides automatic interpretation using the Landis & Koch scale. Users can either launch an interactive graphical interface or use command-line functions for scripting. The package is freely available on CRAN.
0
0
cs.DC 2026-06-05

GPU Newton iteration divides integers at 2^15-2^18 bits near model optimum

by Martin B. Marchioro, Aske N. Raahauge +3 more

On GPU Implementation for Multi-Precision Integer Division

Implementation using shifted-inverse iteration and prefix sums reaches performance close to a multiplication-only cost model on CUDA hardwar

Figure from the paper full image
abstract click to expand
This paper presents the issues arising in implementing a fast integer division algorithm on general purpose GPUs. The algorithm uses a Newton iteration based on the shifted inverse operation, keeping all arithmetic in the integer domain and relying on data-parallel operators. The principal contribution is an efficient GPU/CUDA implementation for integer precisions from $2^{15}$ to $2^{18}$ -- sizes not supported by \cgbn{} division. We propose algorithmic refinements, define a cost model in terms of multiplications, build on prefix sums and previous work on multi-precision multiplication, and present an evaluation showing near-optimal performance relative to the model for the target precision.
0
0
cs.CG 2026-06-05

RedZeD reduces persistent homology to zero differentials for faster pairing

by Chris Kapulkin, Nathan Kershaw

RedZeD: Computing persistent homology by Reduction to Zero Differentials

Active enumeration on Vietoris-Rips filtrations cuts time and memory versus standard algorithms.

Figure from the paper full image
abstract click to expand
We introduce a new algorithm for computing persistent homology of Vietoris--Rips filtrations, which in many cases offers a considerable improvement both in terms of time and memory over the existing implementations of the persistence pairing algorithm. The key innovation, called active enumeration, is made possible by a new theoretical framework of Reduction to Zero Differentials (hence RedZeD) in which to view persistent homology.
0
0
cs.PL 2026-06-04

Low-bit tagging remains fastest for symbolic workloads

by Stephen M. Watt

Look Before You Leap: Checking In on Type Tag Checking

Microbenchmarks on AArch64 and x86-64 show local bit operations beat heap reads for tags while NaN-boxing saves allocation for floats.

Figure from the paper full image
abstract click to expand
Tagging of generic dynamic values is important in symbolic-computation and dynamic-language systems, but the trade-offs change as machine architectures and workloads evolve. In particular, old folklore about boxed values, immediate values, and type tags must be recalibrated from time to time. We revisit the performance of badged object headers, low-bit tagging, and two NaN-boxing layouts on a range of platforms in use today, including AArch64 and x86-64 architectures from different manufacturers. The experiments isolate two distinct effects: the cost avoided by not heap-allocating common scalar values, and the cost avoided by obtaining tag information from the value word rather than by performing a heap read. The results show that several local bit operations are often cheaper than opening a heap object to obtain a tag or small value. Low-bit tagging remains the simplest and usually fastest choice for mostly symbolic workloads, while NaN-boxing is close in access cost and avoids the time and space of heap allocation for ordinary floating-point values.
0
0
cs.AR 2026-06-04

Golden ratio rule fixes exponent widths for all GoldenFloat formats

by Dmitrii Vasilev

GoldenFloat: A Phi-Derived Static-Split Floating-Point Family from GF4 to GF1024 with a Lucas-Exact Integer Identity

The formula matches nine realized formats exactly and generates consistent splits up to 1024 bits without further choices.

abstract click to expand
We present a hardware-oriented description of GoldenFloat (GF), a static-split floating-point family generated by a single closed rule, and three concrete artefacts: (i) an open multi-width RTL generator covering GF4-GF256 with a continuous-integration differential sweep against a correctly-rounded reference; (ii) an integer-backed Lucas-exact accumulator path verified at 500-digit precision for n = 1, ..., 256; and (iii) a GF16 FPGA codec passing a 35-of-35 testbench at 323 MHz on Artix-7 (Xilinx XC7A35T). A format-conformance oracle (Corona) ships in the same repository and is used as the blackbox check in our continuous-integration audit. The rule and its scope. For each total width N >= 4, the exponent width is e = round((N-1)/phi^2) with fraction f = N-1-e and phi = (1+sqrt(5))/2. The rule reproduces the realised exponent widths of nine formats GF4, GF8, GF12, GF16, GF20, GF24, GF32, GF64, GF256 (9/9) and extends consistently to GF128, GF512, GF1024. The rule is positioned alongside posit (2022 Posit Standard), takum (Hunhold 2024, 2025), OCP-MX (Rouhani et al. 2023), and the IEEE P3109 multi-width float draft, all of which are width-spanning families under a parameterised rule. We make no per-rung accuracy or superiority claim against any of them. What is open. The breadth/toolchain-coherence framing is recorded as an open conjecture with a pre-registered falsification path: a matched-substrate FPGA experiment and a matched-budget software ablation. A falsification ledger (FL-002) records the open questions and the experiments that would settle them. An RTL-correctness erratum dated 2026-05-31 is reported in Section 5.5; the fabricated TTSKY26b dies carry the defective multiplier portfolio, and the corrected generator is the regeneration baseline.
1 0
0
math.NA 2026-06-04

Tight bounds give optimal Lipschitz fits with monotonicity

by Gleb Beliakov

Fitting scattered data with optional monotonicity constraints on GPU: LipFit package

Upper and lower envelopes on scattered points produce continuous approximations that obey ordering rules and run on graphics hardware.

Figure from the paper full image
abstract click to expand
This paper presents a method of multivariate scattered data interpolation and approximation that produces optimal Lipschitz-continuous approximation, subject to the desired monotonicity constraints. This method relies on tight upper and lower approximations to the data, and is similar in its spirit to the nearest-neighbour approximation but does not suffer from discontinuities. Local Lipschitz interpolation and Lipschitz smoothing are also presented. This approach falls under the umbrella of instance-based approximation with no training phase, and it is suitable for GPU-based parallelisation. A Python GPU-friendly package LipFit which implements the methods discussed is discussed.
0
0
cs.MS 2026-05-29

FP8 and Kulisch enable FP64 FFT near memory limit on B300

by Satoshi Matsuoka

FP8 is All You Need (Part 2): Efficient Ozaki-Bailey Style FFT Through Tensor-core Garner Reformulation and Kulisch Escape Route

Ozaki-Bailey projects 18 ms for 1024 cubed at full precision versus 12.9 ms memory roof

Figure from the paper full image
abstract click to expand
NVIDIA's Blackwell Ultra (B300) cuts FP64 vector throughput to ~1.3 TFLOPS per GPU, roughly 30x below B200 and well below the level at which bandwidth-limited FP64 workloads stay memory-bound. The Ozaki Scheme II framework recovers FP64-equivalent throughput by routing dense matrix multiply through FP8 tensor cores with a mantissa-sliced Chinese-remainder reconstruction. A companion Part (1) paper covers dense GEMM, batched GEMV, stencils, and SpMV; this paper adds the fifth canonical primitive, the 3-D FFT. We present Ozaki-Bailey FFT, an emulated 3-D FFT via the Bailey six-step decomposition with both 1-D FFT GEMMs on FP8 tensor cores. Bailey's small inner factor k ~ sqrt(N) (k=32 for N=1024) puts the kernel in the regime k << r^2, where the third TME parameter gamma (reconstruction latency) binds rather than amortising. Garner reconstruction splits into Phase A (inner products on FP8/INT8 tensor cores, ~1 ms for 1024^3 on B300) and Phase B (per-output reduction). We identify Kulisch fixed-point complete arithmetic as a Phase B reformulation that keeps full FP64 accuracy while running entirely on the INT32 SIMT pipe. We derive closed-form bandwidth-parity floors. The native FP64 floor is 1.56*B_HBM (12.5 TF at 8 TB/s): B300's 1.3 TF sits ~10x below, Rubin's 33 TF within 4%. The Kulisch escape route needs an INT32 sub-floor 8.25*B_HBM and an FP8 floor 170*B_HBM; B300 meets both. The projection is ~18 ms for 1024^3 at full FP64, essentially the 12.9 ms memory roof. A GPU meets memory-roof FFT parity if it satisfies either the native floor or both Kulisch floors. If the projection holds in practice, B300 becomes viable for full-FP64 FFT through software alone, motivating a libKulisch library and benchmark campaign.
0
0
cs.MS 2026-05-29

Library supplies exact MLE for sixteen HMM emission types

by Gary Wolfman

libhmm: A Modern C++20 Library for Hidden Markov Models with Correct MLE Emission M-Steps

Replaces method-of-moments approximations in Baum-Welch with ECME and Newton-Raphson for Gamma, Beta, Weibull and others.

Figure from the paper full image
abstract click to expand
We describe libhmm, a C++20 library for Hidden Markov Model parameter estimation, sequence decoding, and model selection. libhmm addresses two gaps in existing software: the absence of a well-maintained, zero-dependency C++ HMM library suitable for embedding in production systems, and the widespread use of method-of-moments (MOM) approximations in the emission distribution M-step of the Baum-Welch algorithm. The library implements correct maximum likelihood estimators for sixteen continuous and discrete emission distributions, including an ECME algorithm for the location-scale Student-t distribution, Newton-Raphson maximization for Gamma, Beta, Weibull, and Negative Binomial distributions, and the von Mises distribution for circular data. All forward-backward and Viterbi calculations operate in full log-space. SIMD acceleration is provided for AVX-512, AVX2, SSE2, and ARM NEON via compile-time dispatch with scalar fallback. Python bindings are available via the companion package pylibhmm. We compare libhmm against established C and C++ HMM libraries and against published R reference packages on five real-data benchmarks, and discuss the architectural tradeoffs made in the design.
0
0
math.NA 2026-05-25

Statistical solution converges for Mach 2000 jet at rate 0.5

by Stephan Simonis, Gauthier Wissocq

Computing statistical solutions of a Mach 2000 astrophysical jet

Ensemble of 1000 simulations shows probability measure stabilizing despite chaotic individual flows diverging

Figure from the paper full image
abstract click to expand
The simulation of extreme Mach astrophysical flows is traditionally viewed through the lens of deterministic positivity-preserving schemes. However, due to Kelvin--Helmholtz instabilities and shock anomalies, the multi-dimensional Euler equations admit a variety of non-unique entropy solutions in turbulent regimes. Here, we computationally explore the limits of weak-strong uniqueness of a Mach 2000 jet by defining the statistical solution as the pushforward of a probability measure through a vectorial lattice Boltzmann method operator. Utilizing optimized CUDA kernels, we compute an ensemble of 1000 Monte Carlo samples across a sequence of highly refined spatial grids of up to 3.2 million cells and subsequently post-process the empirical measures via memory-mapped CPU streaming. We contrast the strong sample-wise $L^1$ error divergence with the convergence of the probability measure in the 1-point Wasserstein distance via empirical Cauchy rates. Our results demonstrate that while individual flow realizations physically diverge due to chaotic shear-layer instabilities, the statistical solution converges to an admissible limit measure at a rate of 0.5. Consequently, we provide numerical evidence that the statistical solution to the considered problem is non-Dirac and remains stable in the extreme compressible regime.
0
0
quant-ph 2026-05-25

Julia package computes exact Haar integrals symbolically

by {L}ukasz Pawela, Zbigniew Pucha{l}a

IntegrateUnitary.jl: A Julia package for symbolic integration over Haar measures

Handles Weingarten calculus for unitary, orthogonal, symplectic and other groups with symbolic dimension support.

Figure from the paper full image
abstract click to expand
Symbolic integration over the Haar measure of compact groups is a computational cornerstone in quantum information science and random matrix theory. We present \texttt{IntegrateUnitary.jl}, a comprehensive Julia package for computing exact expectations of polynomial functions over a wide range of compact groups ($U(d)$, $O(d)$, $Sp(d)$, and $SU(d)$ for balanced polynomials), circular and Gaussian ensembles, Ginibre ensembles, permutation groups, random pure states, and unitary $t$-designs. The package provides a fully open-source implementation of the Weingarten calculus and Wick contractions with broad symbolic-$d$ support for entry-wise and trace-polynomial integrals, while selected workflows currently require concrete integer dimensions (including higher pure trace moments $|\mathrm{tr}(U)|^{2k}$ for $k > 1$ and HCIZ with \texttt{SymbolicMatrix} inputs, and direct matrix-valued integration of \texttt{SymbolicMatrix}/\texttt{SymbolicMatrixProduct} expressions), automatic asymptotic expansions, a high-level symbolic trace interface that reconstructs Weingarten graphs from index-free expressions, and a bridge to \texttt{ITensors.jl} for tensor network averaging. We discuss the underlying algorithms, including the Murnaghan-Nakayama rule and symplectic-orthogonal duality, and demonstrate that the package efficiently handles high-degree moments and quantum information metrics.
0
0
math.CO 2026-05-22

Negative counts speed Ehrhart polynomials for Gelfand-Tsetlin polytopes

by Per Alexandersson

Fast computation of Ehrhart polynomials of Gelfand--Tsetlin polytopes via Macdonald reciprocity

Ehrhart-Macdonald reciprocity turns many negative evaluations into quick or zero strict-pattern counts, and adaptive choice yields practical

abstract click to expand
We describe an efficient method for computing the Ehrhart polynomial of Gelfand--Tsetlin polytopes arising from Kostka coefficients. The key idea is to exploit Ehrhart--Macdonald reciprocity: evaluating the Ehrhart polynomial at negative integers reduces to counting \emph{strict} Gelfand--Tsetlin patterns, which are often zero or very small for low dilations. Combined with an adaptive strategy that chooses the cheapest evaluation point (positive or negative) at each step, this yields substantial practical speedups compared to general-purpose polytope software. We benchmark against $\mathtt{OSCAR}$/$\mathtt{polymake}$, and illustrate the broader applicability of the method through order polytopes and permutation posets. The implementation is available in the Rust \texttt{kostka} package, with related optimizations also incorporated in the new \texttt{lrcalc-rs} replacement for \texttt{lrcalc}.
0
0
cs.DC 2026-05-22 2 theorems

Framework computes determinants securely on edge servers

by Prajwal Panth

Secure and Parallel Determinant Computation for Large-Scale Matrices in Edge Environments

Distortion method hides matrix data from untrusted servers yet returns the exact original determinant after parallel processing.

Figure from the paper full image
abstract click to expand
The advent of edge computing has enabled resource-constrained clients to delegate intensive computational tasks to distributed edge servers, especially within Internet of Things (IoT) environments. Among such tasks, Matrix Determinant Computation (MDC) remains critical for applications in control systems, cryptography, and machine learning. However, the cubic complexity of traditional determinant algorithms makes them unsuitable for real-time processing in constrained edge scenarios. We propose a Secure Parallel Determinant Computation (SPDC) framework, which provides strong security guaranties, including privacy-preserving MDC, across N distributed edge servers. The framework achieves privacy through Composite Element Distortion (CED) - a lightweight encryption method that combines Element-wise Obfuscation (EWO) and the Panth Rotation Theorem (PRT) to conceal both structural and numerical matrix content while preserving determinant properties. Parallel LU decomposition is used to distribute encrypted matrix blocks across an arbitrary number of untrusted edge servers, enabling efficient and scalable determinant computation. A one-way communication model further reduces coordination overhead by eliminating inter-server interactions. To ensure result integrity with minimal client burden, we further introduce two verification algorithms: Q_2, a probabilistic scalar method, and Q_3, a deterministic and low-complexity alternative. Mathematical analysis demonstrates that the proposed framework provides strong privacy and security guaranties, low computational overhead, and deployment flexibility - making it well-suited for secure, scalable, and real-time MDC in distributed edge-assisted systems.
0
0
cs.MS 2026-05-21 2 theorems

One toolbox solves both polynomial systems and multiparameter eigenvalue problems

by Christof Vermeersch, Bart De Moor

Solving Multivariate Polynomial Systems and Rectangular Multiparameter Eigenvalue Problems with MacaulayLab

A shared numerical linear algebra method works regardless of basis or monomial order and handles infinite positive-dimensional solution sets

Figure from the paper full image
abstract click to expand
We present the Matlab toolbox MacaulayLab, which implements numerical linear algebra algorithms for solving multivariate polynomial systems and rectangular multiparameter eigenvalue problems. Its structure and functionality are the result of several years of research and algorithmic development. We demonstrate how the software works and compare its performance with other software packages, such as PNLA, PHCpack, and MultiParEig. Some core features of MacaulayLab are the fact that it solves two key problems via one common approach, works independently of the chosen polynomial basis and monomial order, and is capable of dealing with positive-dimensional solution sets at infinity. The toolbox (including its future updates) and a large collection of test problems are freely available online.
0
0
cs.MS 2026-05-20 1 theorem

Smith normal form computation scales with matrix multiplication

by Ziwen Wang, Stavros Birmpilis +2 more

A C implementation of the Smith massager algorithm

Practical C code for the 2020 Las Vegas algorithm runs in time of one BLAS multiply for n up to 10007

Figure from the paper full image
abstract click to expand
We describe a C implementation of the Las Vegas algorithm of Birmpilis, Labahn and Storjohann from 2020 for computing the Smith normal form of a nonsingular integer matrix. The algorithm computes a Smith massager for the input matrix using $O(n^{\omega}\, \B(\log n + \log \|A\|)\, (\log n)^2)$ bit operations, which is softly equivalent to the cost of multiplying two matrices of the same dimension and entry size. We describe the key implementation techniques that bridge the gap between the theoretical algorithm and practical performance, including BLAS-accelerated modular arithmetic via the Residue Number System and an adaptive batching scheme that collapses the theoretical $O(\log n)$ iterations to $O(1)$ in practice. Experiments on matrices of dimension up to $n = 10007$ show that the implementation's running time scales proportionally to that of a single BLAS matrix multiplication, with both exhibiting the same effective growth rate on a log-log plot.
0
0
cs.MS 2026-05-19 1 theorem

Python package detects modes in distributions 3-10 times faster than R

by Ruiyu Zhang, Qihao Wang

critband: A Python Package for Critical Bandwidth Analysis of Multimodal Distributions

critband applies critical bandwidth search with FFT acceleration to give stable mode counts for separated cases across scientific data types

abstract click to expand
Multimodal density estimation is a fundamental problem in scientific computing. Determining the number of modes in a distribution is a core numerical challenge with applications across ecology, economics, genomics, and astronomy. While the R ecosystem provides mature tools through the multimode package, the Python ecosystem has lacked an equivalent cohesive implementation. We present critband, a Python package for critical bandwidth bimodality detection based on Silverman's kernel density approach. The package implements critical bandwidth search with a robust bracketed mode-count solver and FFT-accelerated KDE, and provides additional features including k-mode detection, component decomposition, bimodality strength quantification, and excess mass estimation. Validation against twelve benchmark cases spanning separation regimes, unequal variances, unequal weights, and small sample sizes shows stable estimates for clearly separated cases and expected instability for boundary cases. Performance benchmarks show critband is typically 3-10 times faster per case than R's modetest() in the tested setup.
0
0
cs.PL 2026-05-19 2 theorems

Egglog turns higher-order math models into natural LaTeX and fast constraints

by Hiromi Ishii

Optimizing Optimizations, Declaratively: Optimizing the Higher-Order Functions in Mathematical Optimization with egglog

Equality saturation with datalog rules and Henkin constants reduces domain-set detection from minutes to seconds while restoring readablecom

Figure from the paper full image
abstract click to expand
We present two applications of egglog to mathematical optimization in JijModeling 2, a mathematical modeller whose internal representation is based on simply typed $\lambda$-calculus. First, we use egglog to improve $\LaTeX$ output for mathematical models expressed with higher-order functions. Python comprehensions are desugared into stream operations such as $\textsf{map}$, $\textsf{flat_map}$, and $\textsf{filter}$; emitting these terms directly produces unnatural mathematical notation. We reconstruct comprehension syntax by \emph{ensugaring} higher-order terms and use equality saturation with a custom cost model to minimize temporary variable rebindings. Second, we use egglog as a declarative engine for \emph{constraint detection}, extending the previous egg-based approach presented at EGRAPHS '25. Egglog's datalog-style rules let us express multi-step detection logic directly, without external Rust orchestration code. We encode parametrized constraints using \emph{Henkin-like constants} and propagate side conditions on subterms and indices through egglog facts. Finally, we show that the same ensugaring procedure also reduces large domain-set conditions before saturation, turning a problematic detection case from minutes or nontermination into a few seconds. Through these topics, we want to provide an example of an industrial application of egglog, demonstrate the trick to propagate the constraints using the ideas from mathematical logic, and show the importance of optimizing \emph{premises} of egglog rules to get practical performance in egglog programs.
0
0
cs.CE 2026-05-18 Recognition

Hybrid optimizer improves spatial packaging by over 10 percent

by S. Westerhof, T. Hofman

A Hybrid Optimization Framework for Spatial Packaging of Interconnected Systems

Method using ball decomposition and gradient refinement matches analytical optima within 2 percent on benchmarks.

Figure from the paper full image
abstract click to expand
This paper presents an optimization framework for Spatial Packaging of Interconnected Systems with Physical Interactions (SPI2) that addresses the geometric challenges of three-dimensional component placement and routing. While SPI2 generally includes physical interactions, this study isolates the spatial optimization aspect to evaluate placement and routing performance independently. The framework integrates the Maximal Disjoint Ball Decomposition (MDBD) for geometric abstraction with a hybrid optimization strategy that combines stochastic initialization and gradient-based refinement with interior point optimization. It is formulated to handle the nonlinear, non-convex, and continuous characteristics of spatially coupled design problems. The proposed framework is evaluated against a use case from prior SPI2 research and tested with a newly introduced benchmark that enables verifiable assessment of optimization performance. Results indicate that the presented method achieves more than a 10% improvement over existing SPI2 implementations and converges to spatially analytical optima across various benchmark scenarios. Benchmark experiments show solution accuracy of 0.6-2% relative to the ground truth.
1 0
0
cs.CE 2026-05-18 2 theorems

Optimizer packs components into non-convex custom shapes

by S. Westerhof, T. Hofman

Spatial Optimization of Interconnected Systems in Non-Convex Design Spaces

Differentiable boundary evaluation enables simultaneous routing and inertia optimization while preserving geometric feasibility.

Figure from the paper full image
abstract click to expand
This paper presents a spatial optimization methodology that extends the Spatial Packaging of Interconnected Systems with Physical Interaction (SPI2) framework to support arbitrary, non-convex design boundaries. We introduce a smooth, differentiable inside-outside evaluation for components represented using the Maximal Disjoint Ball Decomposition (MDBD) method. The framework also incorporates center-of-gravity and moment-of-inertia calculations directly into the optimization, and provides an end-to-end computer-aided design (CAD) workflow for importing components and reconstructing the optimized assembly. The method is demonstrated on a fictional aircraft auxiliary unit. Results show that the optimizer can place multiple interconnected components within a custom geometry while simultaneously handling routing and physics-based objectives. The approach maintains geometric feasibility within numerical tolerance and illustrates the potential of MDBD-based SPI2 methods for practical engineering design applications.
1 0
0
cs.LO 2026-05-18 2 theorems

Lean formally verifies full BET surface area workflow

by Ejike D. Ugwuanyi, Colin T. Jones +2 more

LeanBET: Formally-verified surface area calculations in Lean

The pipeline proves regression coefficients minimize error, window enumeration is sound, and results match reference to machine precision.

Figure from the paper full image
abstract click to expand
The Brunauer--Emmett--Teller (BET) method is a standard tool for estimating surface areas from adsorption isotherms, yet practical implementations involve multiple algorithmic steps whose correctness is rarely made explicit. In this work, we present a fully executable and formally verified BET analysis pipeline implemented in the Lean~4 theorem prover. Our formalization covers the complete BET Surface Identification (BETSI)-style workflow, including window enumeration, monotonicity checks, knee selection, and linear regression. We carry out computations in floating-point arithmetic and develop the corresponding correctness proofs over the real numbers, using a shared polymorphic implementation that supports both. On the proof side, we show that the regression coefficients returned by the algorithm agree with their specification-level definitions and minimize the least-squares error under the stated assumptions. We also formalize the algebraic derivation of the BET linearized expression and connect that result directly to the executable analysis pipeline. We further prove that the window enumeration is sound and complete, and that the admissibility checks and knee-based selection satisfy their formal specifications. We evaluate the implementation against the BETSI reference method on benchmark adsorption isotherms. Compared to BETSI, LeanBET agrees to machine precision for 18 of the 19 isotherms, with only a 0.03\% deviation for the UiO-66 dataset. This demonstrates that a scientific computing workflow can be built in Lean, yielding both formal verification guarantees and numerical agreement with an established Python reference implementation.
0
0
cs.DC 2026-05-18 Recognition

Parallel code speeds star-M SVD compression for big datasets

by Md Taufique Hussain, Grey Ballard +4 more

High-Performance Star-M SVD for Big Data Compression

New implementation delivers optimal compression of large scientific data with matrix-like guarantees.

Figure from the paper full image
abstract click to expand
In the era of big data, effectively compressing large datasets while performing complex mathematical operations is crucial. Tensor-based decomposition methods have shown superior compression capabilities with minimal loss of accuracy compared to traditional matrix methods. Under the star-M tensor framework, tensors can be decomposed in a matrix-mimetic way, including using the star-M SVD. This tensor SVD has optimality guarantees and has shown exceptional performance on specific types of data, but software implementations have been mostly limited to productivity-oriented languages. In this work, we present our development of a shared-memory parallel, high-performance solution designed to efficiently implement the underlying algorithms. This software will enable optimal compression of extensive scientific datasets, paving the way for enhanced data analysis and insights.
0
0
cs.MS 2026-05-18 Recognition

SIMD algorithms enable first correctly rounded vector math library

by Cristina Anderson, Marius Cornea +2 more

Correctly Rounded Functions For Vector Applications: A Performance Study

Single-precision designs integrated into CPU library form core of release planned for mid-2026, with GPU versions for cross-platform results

Figure from the paper full image
abstract click to expand
Following recent interest in correctly rounded math library functions (as currently recommended by the IEEE 754 standard), we have designed several SIMD algorithms for one-input single precision functions and integrated them into our CPU math library; these will form the core of the first correctly rounded vector math library, to be available to users in mid-2026. To take advantage of the cross-platform bitwise reproducibility afforded by correct rounding, we adapted and evaluated a few SIMD implementations on graphics processing units (GPU). In addition, we designed and evaluated proof-of-concept SIMD implementations of two correctly rounded double precision functions.
0
0
cs.MS 2026-05-14

Interior point methods run on GPUs by compressing sparse problems

by Slaven Peles, Kalyan S. Perumalla +5 more

Porting the Nonlinear Optimization Library HiOp to Accelerator-Based Hardware Architectures

The HiOp library reformulates each iteration so the entire solver fits on accelerator hardware for optimal power flow cases.

Figure from the paper full image
abstract click to expand
While interior point methods have been the centerpiece of nonlinear programming tools used in science and engineering, their reliance on linear solvers that can tackle sparse symmetric indefinite and highly ill-conditioned problems made it difficult to implement them effectively on hardware accelerators. At this time, there are few sparse linear solvers that can be used in this context. Here, we present a novel formulation of an interior point method implemented in our HiOp library, which is designed to be able to run entirely on hardware accelerators. This formulation avoids dependence on sparse solvers altogether, which is achieved by compressing the underlying sparse linear problem into a dense one of manageable size. We demonstrate feasibility of this approach and provide a baseline for future interior point method implementations on hardware accelerators. Our investigation is motivated by problems arising in optimal power flow analysis in power systems engineering and our approach is tailored to the broad class of problems arising in that important domain. We also demonstrate utility of modern programming models based on performance portability libraries, namely, Umpire and RAJA. We discuss trade-offs between performance, portability and development cost in the solution space for this non-linear optimization problem. As a result of this research, we demonstrate for the first time that interior point methods for sparse problems can be efficiently realized on modern computing systems where more than 90% of processing power is in GPUs.
0
0
stat.CO 2026-05-14

Python library joins simulators with time-average diagnostics

by Ihor Kendiukhov

Ergodicity Library: A Python Toolkit for Stochastic-Process Simulation, Time-Average Diagnostics, and Agent-Based Experiments

Reduces custom code needed for non-ergodic dynamics and uncertainty experiments

Figure from the paper full image
abstract click to expand
ergodicity is an open-source Python library for computational work on stochastic dynamics, with particular emphasis on non-ergodicity, time-average behavior, heavy-tailed processes, and decision making under uncertainty. The package brings together three layers that are often split across ad hoc scripts: process definitions and simulators, analysis and fitting tools, and agent-based experimentation. This article documents the implemented software rather than presenting new stochastic theory. We describe the package architecture, the supported process families, the analysis workflow, and the practical boundaries of the current implementation. We also provide fully reproducible examples covering heavy-tailed ensemble spread, multiplicative Levy growth diagnostics, adaptive memory mean reversion, preasymptotic fluctuation analysis, and partial stochastic differential equation simulation. The package is positioned as an integration layer on top of the scientific Python stack, reducing the amount of glue code required to move from process specification to diagnostics and comparative experiments.
0
0
eess.SY 2026-05-13 2 theorems

Docker container makes Basilisk GN&C simulations reproducible

by Anubhav Gupta

Basilisk and Docker for Reproducible GN&C Simulation: A Workflow Reference

The workflow packages the full environment so identical spacecraft control tests run on any machine without dependency conflicts.

Figure from the paper full image
abstract click to expand
Basilisk is an open-source astrodynamics simulation framework widely used for spacecraft guidance, navigation, and control (GN&C) research and development. Despite its flexibility and computational capabilities, configuring Basilisk consistently across heterogeneous development environments presents practical challenges due to dependency management, operating system compatibility, and software configuration requirements. This paper presents a Docker-based containerization workflow for Basilisk that encapsulates the complete build environment, dependencies, and simulation infrastructure within a portable container image. The workflow is demonstrated through a progression of simulation scenarios of increasing complexity, from standalone orbital dynamics scripts to BSKSim-based attitude dynamics and control simulations with Monte Carlo analysis. The BSKSim class hierarchy, dynamics model architecture, flight software implementation, and scenario execution patterns are described in detail. The presented workflow provides a self-contained implementation reference for GN&C engineers and researchers seeking reproducible and portable Basilisk simulation environments. This work expands upon a workshop presentation delivered at the 46th Rocky Mountain AAS GN&C Conference, February 2024, available at https://doi.org/10.5281/zenodo.15008785.
0
0
quant-ph 2026-05-13 1 theorem

Compiler partitions multi-QPU circuits to cut port overflow and link congestion

by Soumyadip Sarkar, Subhasree Bhattacharjee

QuPort: Topology-, Port-, and Congestion-Aware Compilation for Modular Multi-QPU Quantum Systems

TPCCAP balances weighted cut distance, port usage, and routed loads on a three-level topology model so cross-device traffic does not create瓶

Figure from the paper full image
abstract click to expand
Modular quantum processors require a compiler to reason about two resources at the same time: local device connectivity and communication across QPUs. A mapping that is acceptable on a single coupling graph may be unsuitable for a modular machine if it creates excessive cross-QPU traffic, concentrates that traffic on a small number of interconnect links, or assigns many boundary qubits to a QPU with few communication ports. This paper presents QuPort, a Python and Qiskit-based compilation framework that studies this setting through an explicit three-level model: a weighted logical interaction graph, a directed physical coupling map, and an undirected QPU-level interconnect graph. The main partitioning method, TPCCAP, optimizes the implemented objective formed by weighted cut distance, communication-port overflow, and routed link-load congestion. The framework also includes heavy-edge clustering, balanced greedy partitioning, simulated-annealing refinement, communication-port-aware layout, extraction of remote two-qubit operations, local-only routing of per-QPU circuits, and topology-aware schedule estimation. The model is a compiler-level abstraction. It does not claim a calibrated hardware runtime or an implementation of a physical remote-gate protocol.
0
0
cs.CE 2026-05-12 Recognition

Distributed NUFFT scales to 1024 GPUs on mixed hardware

by Paul Fischill, Andreas Adelmann +1 more

A Performance-Portable, Massively Parallel Distributed Nonuniform FFT

Kokkos code matches specialized single-GPU speed while enabling billion-particle spectral plasma runs on three supercomputers

Figure from the paper full image
abstract click to expand
The nonuniform fast Fourier transform (NUFFT) enables spectral methods for problems with irregularly spaced samples, with applications in medical imaging, molecular dynamics, and kinetic plasma simulations. Existing implementations are limited to shared-memory execution, restricting problem sizes to what fits on a single node. We present the first distributed, performance-portable NUFFT for heterogeneous supercomputers. Our Kokkos-based implementation runs without modification on NVIDIA and AMD GPUs. We develop multiple spreading and interpolation kernels optimized for different accuracy requirements and architectures. Our spreading kernels match or exceed the single-GPU throughput of the state-of-the-art CUDA-based NUFFT library cuFINUFFT at production particle densities, while our Kokkos-based implementation additionally supports AMD GPUs. Strong scaling experiments on Alps (NVIDIA GH200), JUWELS Booster (NVIDIA A100), and LUMI (AMD MI250X) demonstrate scaling up to 1024 GPUs. At scale, the distributed FFT is a significant part of the total runtime, making higher NUFFT accuracy less expensive. We apply the method to massively parallel Particle-in-Fourier simulations of Landau damping with up to $1024^3$ Fourier modes and 8.6 billion particles on Alps, JUWELS, and LUMI, demonstrating that distributed NUFFTs enable kinetic plasma simulations at resolutions previously inaccessible to spectral particle methods.
0
0
math.OC 2026-05-12 2 theorems

Riemannian L-BFGS handles Euclidean bounds on manifolds

by Mateusz Baran, Ronny Bergmann +1 more

A Riemannian quasi-Newton algorithm for optimization with Euclidean bounds

Tangent-space updates plus adapted Cauchy point strategy match classical performance on Euclidean problems and deliver large speedups on two

Figure from the paper full image
abstract click to expand
We propose a Riemannian limited-memory BFGS method for optimization problems with Euclidean bounds. The method combines a limited-memory quasi-Newton update in the tangent space with a Riemannian adaptation of the generalized Cauchy point strategy from classical L-BFGS-B, enabling efficient handling of Euclidean bounds while exploiting the geometric structure of the optimization domain. This setting is important in several applications, including covariance matrix estimation with bounded variance, neuroimaging, EEG signal classification, and other signal processing or computer-vision tasks that couple manifold variables with constrained Euclidean parameters. We provide a generic algorithmic framework and an implementation of the algorithm in the Manopt.jl library. Numerical experiments on benchmark problems indicate only minor reduction in performance on Euclidean problems compared to the classical L-BFGS-B method, while outperforming interior-point methods. Furthermore, the algorithm was tested on two mixed manifold and bounded Euclidean problems: amplitude-limited blind source separation with Gaussianity penalization and bounded-variance maximum likelihood common principal components analysis. The proposed method outperforms existing methods by several orders of magnitude.
1 0
0
cs.MS 2026-05-11 Recognition

GPU solver speeds up entropic optimal transport calculations

by Yixuan Qiu

cuRegOT: A GPU-Accelerated Solver for Entropic-Regularized Optimal Transport

cuRegOT's optimizations deliver major performance gains over existing GPU methods while ensuring algorithmic convergence.

Figure from the paper full image
abstract click to expand
Optimal transport (OT) has emerged as a fundamental tool in modern machine learning, yet its computational cost remains a significant bottleneck for large-scale applications. While harnessing the massive parallelism of modern GPU hardware is critical for efficiency, the de facto standard Sinkhorn algorithm, despite its ease of parallelization, often suffers from slow convergence in challenging problems. More recently, the sparse-plus-low-rank quasi-Newton method offers a balance between convergence rate and per-iteration complexity; however, its efficiency on GPUs is severely hindered by the serial nature of sparse matrix symbolic analysis and irregular memory access patterns. To bridge this gap, we present cuRegOT, a high-performance GPU solver tailored for entropic-regularized OT. We introduce a suite of algorithmic and architectural optimizations, including an amortized symbolic analysis strategy to mitigate CPU bottlenecks, an asynchronous Sinkhorn iterates generation mechanism, and a fused kernel for bandwidth-efficient gradient evaluation. These strategies are backed by rigorous theoretical guarantees ensuring algorithmic convergence. Extensive numerical experiments demonstrate that cuRegOT achieves significant speedups over state-of-the-art GPU-based solvers across a variety of benchmark tasks.
0
0
q-bio.QM 2026-05-11 1 theorem

R package MeTime stores all metabolomics steps in one container

by Bharadwaj Marella, Patrick Weinisch +6 more

MeTime: An R package for reproducible longitudinal metabolomics data analysis

Longitudinal studies keep data, metadata and results together so workflows stay transparent and reports generate automatically.

abstract click to expand
MeTime is an opensource R package for reproducible analysis of longitudinal metabolomics data. It builds upon a central S4 container, metime_analyser, that stores multiple datasets, associated metadata and analysis outputs, enabling unified handling of complex longitudinal studies. Analyses are constructed by piping modular functions, beginning with data transformations (mod_), followed by calculations (calc_), and optional meta-analysis (meta_), so entire workflows remain transparent and easy to modify. MeTime wraps numerous existing methods within a consistent interface, including sample and metabolite distributions, correlation and distance matrices, dimensionality reduction (PCA, UMAP, tSNE), random forest imputation and feature selection via Boruta, eigenmetabolites and WGCNA based clustering, conservation index analysis, regression models (linear, mixed effects, and generalized additive), and partial correlation networks. By retaining all intermediate results and provenance within the container, MeTime facilitates iterative exploration and ensures reproducible reporting via automatically generated HTML and PDF outputs. Comprehensive user guides, case studies and reference documentation accompany the package, making MeTime a versatile platform for longitudinal omics workflows.
0
0
cs.DC 2026-05-08 Recognition

FalconGEMM exceeds GEMM speeds by 7-17% via lower-complexity algorithms

by Honglin Zhu, Jiaping Cao +11 more

FalconGEMM: Surpassing Hardware Peaks with Lower-Complexity Matrix Multiplication

The framework automates LCMA selection and group-parallel execution to beat cuBLAS and AlphaTensor on GPUs and CPUs.

Figure from the paper full image
abstract click to expand
Peak breaking Matrix Multiplication is a promising technique to improve the performance of DL, especially in LLM training and inference. We present FalconGEMM, a cross-platform framework that automates the deployment, optimization, and selection of Lower-Complexity Matrix Multiplication Algorithms (LCMAs) across diverse hardware. There are three key innovations: (1) a Deployment Module that enables portable execution across various hardware and input configurations through code generation; (2) an Execution Module with Group-Parallel Optimizations that maximizes on-chip data reuse, utilizes parallel resources, and reduces bandwidth overhead; and (3) a Decision Module featuring a lightweight analytical performance model to select the optimal strategy based on matrix shapes and hardware profiles. Extensive evaluation is conducted on LLM workloads across GPU (H20, A100) and CPU (ARM, x86) architectures with multiple data types. FalconGEMM succeeds in delivering peak breaking performance and outperforms GEMM libraries (e.g., cuBLAS, CUTLASS, Intel MKL, etc) by 7.59%-17.85% and LCMA competitors like AlphaTensor by 12.41%-55.61%. Our framework makes the theoretical promise of LCMAs practical for production deployment across the heterogeneous landscape of modern hardware.
0
0
cs.LG 2026-05-07

Gradients computed for DAE simulations with events

by Ion Matei, Maksym Zhenirovskyy +1 more

Differentiable Parameter Optimization for DAEs with State-Dependent Events

Two methods return valid derivatives along the forward event path for systems with algebraic constraints and resets.

Figure from the paper full image
abstract click to expand
Differential-algebraic equations (DAEs) with state-dependent events arise in systems whose continuous dynamics are constrained by algebraic equations and interrupted by mode changes, switching logic, impacts, or state reinitializations. Gradient-based parameter learning for such systems is challenging because algebraic variables are implicitly defined, event times depend on the parameters, and reset maps introduce discontinuities. This paper studies differentiable parameter optimization for semi-explicit DAEs with events. We formulate the learning problem as a constrained least-squares problem with DAE dynamics, algebraic constraints, guard equations, and reset maps. We then develop two complementary gradient-computation strategies. The first is an automatic-differentiation-through-simulation method that solves algebraic variables inside the vector field, differentiates the algebraic solve using the implicit function theorem, and handles events through segmented differentiable integration. The second is an explicit discrete-adjoint method that represents the forward simulation as an event-split residual system and computes gradients by solving for the Lagrange multipliers of smooth-segment and event residuals. The formulation clarifies that residual terms in the adjoint method are equality constraints, not heuristic penalties. We compare the two approaches in terms of gradient interpretation, event-time handling, implementation complexity, and local validity. Both methods provide gradients for the event path selected by the forward simulation and are valid under fixed event ordering and transversal guard crossings.
0
0
stat.AP 2026-05-07

Randompack produces identical random sequences across languages and machines

by Kristján Jónasson

Randompack: Cross-Platform Reproducible Random Number Generation and Distribution Sampling

The C library separates engines from distributions to achieve both faster sampling and bit-for-bit compatibility from the same seeds.

Figure from the paper full image
abstract click to expand
A C library for random number generation, Randompack, is presented. The library implements several modern random number generators (engines), including xoshiro256, PCG64, Philox, ranlux++, and sfc64; 14 continuous distributions including uniform, normal, exponential, gamma, beta, and multivariate normal; raw bit streams, bounded integers, permutations, and sampling without replacement. The engine and the distribution layers are separated so any engine can be used with any distribution. Benchmarks show that Randompack is faster overall than competing libraries, with speedup factors ranging from about 1 to 15 depending on engine, distribution, interface, and platform. A distinguishing feature is reproducibility: with the same seeds Randompack gives compatible results across programming languages, computers, CPU architectures, and compilers. The library includes comprehensive support for parallel simulation. It is accompanied by a comprehensive test suite, benchmarking programs, and example programs. Interfaces to Fortran, Python, Julia, and R have been implemented; their benchmark results are included, although their design and implementation are otherwise outside the scope of the article. Unlike other available C libraries with comparable scope, Randompack is permissively licensed under the MIT license, and it is open source and publicly available through GitHub and conda-forge.
0
0
cs.MS 2026-05-07

CombOL delivers unbiased Boltzmann samples via dynamic precision

by Casper Asbj{o}rn Eriksen, Daniel Merkle

CombOL: a Library for Practical Enumeration and Boltzmann Sampling of Combinatorial Classes

A library takes compact class definitions, derives generating functions, and returns statistically exact samples free of floating-point bias

Figure from the paper full image
abstract click to expand
We present CombOL (Combinatorial Objects Library), an open-source library for the enumeration and Boltzmann sampling of combinatorial classes. Classes can be specified by a concise string syntax, and may depend on an arbitrary number of parameters. CombOL automatically derives the associated generating functions, enabling the generation of counting sequences and the compilation of Boltzmann samplers. The library supports exact and approximate-size Boltzmann rejection sampling with automatic parameter tuning to target specific sizes. In addition to implementing established methods, CombOL contributes a novel early-rejection scheme, as well as guaranteed statistical correctness by dynamically increasing the numerical precision, eliminating bias due to floating-point rounding errors. Through the Python interface, sampled structures can be mapped to application-specific objects, enabling direct sampling of domain objects such as graphs, chemical structure representations, or other complex data types. CombOL is available from PyPI as 'combol' (pypi.org/project/combol). The source code is available at gitlab.com/casbjorn/combol.
0
0
cs.MS 2026-05-05

mrdi file format supports distributed algebraic computing

by Antony Della Vecchia

Interprocess Communication of Algebraic Data

Design decisions in the serialization framework allow fine tuning for specific use cases while enabling interprocess communication.

Figure from the paper full image
abstract click to expand
We discuss implementation details of OSCAR's serialization framework, highlighting the design decisions that allow the fine tuning of serialization methods for specific use cases. In particular, we show how the mrdi file format can be used for distributed computing.
0
0
cs.RO 2026-05-05 3 theorems

Hybrid seam gives exact SE(3) Hessians five times faster

by Frank O. Kuehnel

Exact Higher-Order Derivatives for SE(3) via Analytical/AD Methods

Placing analytical Jacobians up to y=Tx lets the same gradient code produce precise higher-order tensors without finite-difference tuning.

abstract click to expand
Fast prototyping of new SE(3) estimation objectives remains awkward in practice. Modern Lie-group frameworks -- GTSAM, manif, Sophus, SymForce, Ceres -- target first-order workloads through different code-generation and automatic-differentiation strategies, each optimized for a particular seam between hand-derived geometry and generic differentiation. The remaining gap is a compact, AD-safe path from these first-order primitives to exact Hessians, observed-information matrices, and higher-order derivative tensors: the quantities needed for exact Newton steps, observed-information covariance estimates, and covariance correction. This paper presents a hybrid analytical/AD recipe for SE(3) negative log-likelihoods. The practitioner writes the NLL gradient once, generic over a scalar type, and places the analytical/AD seam at the point-action interface y = Tx. Closed-form Lie-group Jacobians are used up to this interface; AD is applied only beyond it. The same source is then instantiated with ordinary floating-point scalars for gradients, vector-seeded dual numbers for exact Hessians in a single forward-mode pass, and nested dual numbers for higher-order derivative tensors. On a representative 6-DoF, 5-landmark SE(3) NLL, the advocated seeded-Hessian path is approximately 5x faster than finite-differencing the AD gradient on this benchmark while matching a nested-AD oracle to machine precision. The implementation adds roughly 70 lines of analytical-Jacobian code over an AD-only baseline. We also identify and fix a removable singularity in the standard SO(3)/SE(3) scalar basis that would otherwise produce NaNs at the origin under seeded AD, and we audit which Lie-group derivative tensors require this stabilized basis. The result is a practical path from rapidly written SE(3) objectives to exact higher-order derivatives, with predictable runtime and no finite-difference tuning.
0
0
quant-ph 2026-05-04

Multi-objective workflow selects quantum strategies reproducibly

by Soumyadip Sarkar

QBalance: A Reproducible Multi-Objective Workflow for Quantum Compilation, Noise Suppression, and Error-Mitigation Strategy Selection

QBalance uses weighted objectives and non-dominated selection to balance compilation, noise suppression, and mitigation choices for nearterm

Figure from the paper full image
abstract click to expand
Near-term quantum workloads are shaped by coupled compilation and execution choices: qubit layout, routing, basis translation, gate suppression, measurement mitigation, shot budget, and artifact reproducibility. This paper analyzes QBalance, a Python workflow library for dataset-level selection among quantum compilation, noise-suppression, and error-mitigation strategies built on the Qiskit ecosystem. The contribution is formulated as a finite multi-objective strategy-selection problem over circuits, backends, and transformation policies. The manuscript derives the implemented weighted objective, non-dominated selection rule, survival-product error proxy, Bayesian linear candidate-ordering surrogate, and distributional diagnostics. It also positions the system relative to established work on Qiskit pass-manager compilation, SABRE-style routing, randomized compiling, dynamical decoupling, zero-noise extrapolation, matrix-free measurement mitigation, circuit cutting, and Thompson sampling. The analysis shows that QBalance provides a reproducible orchestration and artifact model for quantum workflow studies. It also establishes precise limitations: the current bandit mechanism orders candidates but does not reduce the number of candidate evaluations, the custom layout heuristic is greedy and only partially topology-aware, the implemented ZNE helper is parity-centered, and the cutting integration is a hook rather than a full reconstruction pipeline.
0
0
physics.data-an 2026-05-01

FitED makes nonlinear data fitting reproducible with uncertainty tools

by Mustafa Mahmoud Aboulsaad

FitED: A User-Centric, Extensible Software Environment for Robust Peak-Profile and General Functional Data Fitting

A Python desktop app supports standard and custom models while exporting structured results and propagating parameter uncertainties for expe

Figure from the paper full image
abstract click to expand
Reliable parameter extraction from experimental data is essential for quantitative analysis across spectroscopy, diffraction, photoluminescence, chromatography, microscopy, and time-resolved measurements. However, nonlinear fitting often remains difficult to reproduce, especially when complex models, correlated parameters, uncertain derived quantities, and user-dependent fitting choices are involved. We present FitED, a Python-based desktop application for nonlinear fitting of one-dimensional scientific data that combines an accessible graphical interface with a transparent and flexible numerical backend. FitED supports conventional peak profiles, including Gaussian, Lorentzian, Pseudo-Voigt, and exact area-normalized Voigt functions, as well as arbitrary user-defined analytical models for broader experimental applications. The software integrates local and global-search-assisted optimization strategies, automated model initialization, repeated stability testing, parameter-correlation analysis, and covariance-based propagation of uncertainty for derived quantities. By combining interactive usability with uncertainty-aware analysis and structured export of fitting results, FitED provides a practical platform for reproducible and interpretable fitting of experimental data. The software is intended to support both routine analysis and advanced model evaluation while preserving the parameter-level control required by experimental researchers.
0
0
q-fin.CP 2026-04-30

fast-vollib adds PyTorch and JAX backends to implied-volatility calculations

by Raeid Saqur

Fast-Vollib: A Fast Implied Volatility Library for Pythonwith PyTorch, JAX, and CUDA Fused-Kernel Backends

Library matches py_vollib API while shipping vectorized Halley and Jäckel solvers for batched European options

Figure from the paper full image
abstract click to expand
We present fast-vollib, an open-source Python library that provides high-performance European option pricing, implied volatility (IV) computation, and Greeks under the Black-76, Black-Scholes, and Black-Scholes-Merton models. The library is designed as a drop-in alternative to the de-facto-standard py_vollib and py_vollib_vectorized packages, with pluggable PyTorch and JAX execution backends, a CUDA fused-kernel Triton contribution for batched IV workloads, and a compatibility-first public API. In addition to a vectorized Halley-method IV solver, fast-vollib ships an experimental, fully-vectorized implementation of J\"ackel's "Let's Be Rational" (LBR) algorithm with NumPy/Numba, torch.compile, JAX, and Triton single-pass GPU kernels for batched option chains. This note announces the library and describes its public API surface, with source, documentation, and packaging artifacts available at: GitHub (https://github.com/raeidsaqur/fast-vollib), Docs (https://raeidsaqur.github.io/fast-vollib/), PyPI (https://pypi.org/project/fast-vollib/).
1 0
0
stat.ME 2026-04-29

The paper applies Variational Expectation Maximization to fit nonlinear mixed effects…

by Mohamed Tarek, Pedro Afonso

Fitting Large Nonlinear Mixed Effects Models Using Variational Expectation Maximization

Variational Expectation Maximization scales NLME model fitting to over 15,000 population parameters using flexible variational families and…

abstract click to expand
Nonlinear Mixed Effects models (NLME) models are widely used in pharmacometrics and related fields to analyze hierarchical and longitudinal data. However, as the number of parameters and random effects increases, traditional methods for maximizing the marginal likelihood become computationally expensive. This paper explores the Variational Expectation Maximization (VEM) algorithm, a scalable alternative for fitting NLME models. Originally introduced in the context of probabilistic graphical models and later popularized through variational autoencoders, VEM has not been extensively applied to NLME modeling. By leveraging flexible variational families and reverse-mode automatic differentiation, VEM can efficiently maximize the marginal likelihood, scaling to NLME models with over 15,000 population parameters. This work provides a detailed description of VEM, compares it to other NLME fitting algorithms, and highlights its scalability through computational experiments. Using the Pumas statistical software, we fit two test models: 1) a standard warfarin model, and 2) a DeepNLME Friberg model with 15,410 population parameters and 16 random effects. The warfarin model was fitted to completion to demonstrate the correctness of VEM, while the DeepNLME Friberg model was fitted for a limited number of iterations to measure the time per iteration and demonstrate VEM's scalability.
0
0
quant-ph 2026-04-29

Classical emulator of HHL scales only with qubit count

by Reece Robertson, Ameya Bhave

Extending UNIQuE: Quantum Simulation Speedup for the HHL Algorithm

Emulated version grows exponentially with system size alone, unlike state-vector methods that also track eigenvalue magnitude.

Figure from the paper full image
abstract click to expand
In an extension of the Unconventional Noiseless Intermediate Quantum Emulator, this work introduces a classical emulation of the quantum Harrow-Hassidim-Lloyd algorithm for sampling from the solution space of linear systems. The emulated HHL algorithm scales exponentially with the number of qubits required to represent the linear system, which is an advantage over the state vector simulation of the HHL algorithm, which scales exponentially as a function of both the size of the linear system and the magnitude of its largest (scaled) eigenvalue. We benchmark our emulator by comparing it with the Intel Quantum Simulator and demonstrate a runtime advantage for small linear systems.
0
0
cs.LG 2026-04-28 2 theorems

Asymmetric key-value quantization lowers attention KL error at 4 bits

by Paolo D'Alberto

Statistical Inference and Quality Measures of KV Cache Quantisations Inspired by TurboQuant

Treating keys with simpler quantization than values reduces divergence because softmax amplifies key inner-product variance more than value.

Figure from the paper full image
abstract click to expand
We analyse three KV cache quantization schemes under a fair bit budget: \textbf{KV} (scalar MSE baseline), \textbf{KQV} (WHT + MSE on $K$; WHT + MSE + QJL on $V$), and \textbf{QKQV} (WHT + MSE + QJL on both). Starting from the Beta distribution on the hypersphere, we trace how QJL on $K$ inflates inner product variance by $\pi/2$, which softmax amplifies nonlinearly via Jensen's inequality, and we present statistical inference and information metrics to highlight practical differences. Three empirical findings emerge. (1)~At $n=4$ (the practically dominant budget), KQV wins on every measure -- KL divergence, geometric $K$ error, and 6D distance -- across all distributions and ranks tested. (2)~The K--V asymmetry is unconditional: QKQV is consistently worse than KQV in KL divergence at every budget and distribution. (3)~A budget-dependent crossover exists: QKQV achieves better geometric $K$ reconstruction at $n \in \{2,3,5\}$, KQV at $n \in \{4,6\}$, invariant to rank and tail weight -- an open rate-distortion problem. $\mathrm{KL}(p_{\mathrm{ref}} \| p_{\mathrm{quant}})$, K-only by construction, bridges K direction error to routing corruption and output collapse. We present a sufficient condition when the Jensen mechanism amplifies superlinearly through the softmax. At $n \in \{2,3,5\}$, QKQV wins geometrically because this assumption does not bind. At $n=4$, elevated K error and KL divergence for QKQV strongly suggest the Jensen mechanism is the operative cause of the crossover, providing a new perspective and explanation.
1 0
0
cs.MS 2026-04-27

Compile-time templates fuse GPU kernels to beat PyTorch

by Ryan R. Curtin, Marcus Edel +1 more

Fast GPU Linear Algebra via Compile Time Expression Fusion

Bandicoot matches Armadillo's API while generating optimized kernels without JIT overhead or runtime expression costs.

Figure from the paper full image
abstract click to expand
We describe the Bandicoot GPU linear algebra toolkit, a C++ based library that prioritises ease of use without compromising efficiency. Bandicoot's API is compatible with the popular Armadillo CPU linear algebra library, enabling easy transition for existing CPU-based codebases. Unlike other GPU-focused toolkits, Bandicoot uses template metaprogramming to generate fused GPU kernels directly at compile time, yielding efficient kernels that are often able to saturate memory bandwidth. This removes the need for runtime overhead or JIT infrastructure. Empirical results show that Bandicoot outperforms (sometimes by considerable margins) commonly-used linear algebra toolkits including PyTorch, TensorFlow, and JAX.
0
0
cs.MS 2026-04-24

Hybrid JAX-PETSc code beats pure JAX on large micromechanics simulations

by Alberto Cattaneo, M Keith Ballard +2 more

JetSCI: A Hybrid JAX-PETSc Framework for Scalable Differentiable Simulation

JetSCI keeps JAX differentiation and GPU speed while handing the linear algebra to PETSc for distributed-memory efficiency and accuracy.

Figure from the paper full image
abstract click to expand
The rapid rise of scientific machine learning (SciML) has expanded the role of differentiable modeling, surrogate modeling, and data-driven constitutive laws in large-scale simulation. The JAX framework provides an attractive environment for these workflows through automatically differentiable programs, vectorization, GPU acceleration, and while enabling seamless learning of surrogate models. However, large-scale simulation still relies on mature HPC infrastructure. Libraries, such as PETSc, provide scalable MPI-based parallelism, robust linear and nonlinear solvers, and advanced preconditioning capabilities that remain difficult to reproduce in JAX-only workflows. We present JetSCI, a hybrid JAX-PETSc framework that unifies these complementary strengths. JetSCI uses JAX for GPU-parallel differentiable discretizations and PETSc for robust, scalable solution of the resulting systems on distributed-memory architectures, exposing multilevel parallelism through GPU acceleration within nodes and MPI parallelism across nodes. For finite element discretizations of heterogeneous micromechanics problems, JetSCI outperforms JAX-only implementations in efficiency and accuracy.
0
0
cs.DS 2026-04-24

The paper introduces an algorithm that generates random graphs with prescribed expected…

by Gianlorenzo D'Angelo, Riccardo Michielan

Efficient generation of expected-degree graphs via edge-arrivals

An exact O(n+m) generator for rank-1 expected-degree graphs uses edge arrivals in a multigraph followed by projection to recover the…

abstract click to expand
We study the efficient generation of random graphs with a prescribed expected degree sequence, focusing on rank-1 inhomogeneous models in which vertices are assigned weights and edges are drawn independently with probabilities proportional to the product of endpoint weights. We adopt a temporal viewpoint, adding edges to the graph one at a time up to a fixed time horizon, and allowing for self-loops or duplicate edges in the first stage. Then, the simple projection of the resulting multigraph recovers exactly the simple Norros--Reittu random graph, whose expected degrees match the prescribed targets under mild conditions. Building on this representation, we develop an exact generator based on \textit{edge-arrivals} for expected-degree random graphs with running time $O(n+m)$, where $m$ is the number of generated edges, and hence proportional to the output size. This removes the typical vertex sorting used by widely-used fast generator algorithms based on \textit{edge-skipping} for rank-1 expected-degree models, which leads to a total running time of $O(n \log n + m)$. In addition, our algorithm is simpler than those in the literature, easy to implement, and very flexible, thus opening up to extensions to directed and temporal random graphs, generalization to higher-order structures, and improvements through parallelization.
0
0
cs.DC 2026-04-22

HyperLogLog skips exact counts for faster GPU SpGEMM

by Yifan Li, Giulia Guidi

Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU

Lightweight estimation replaces the symbolic pass to cut runtime 1.4x-2.8x on A100 and H100 across square and rectangular matrices

Figure from the paper full image
abstract click to expand
In computational science and data analytics, many workloads involve irregular and sparse computations that are inherently difficult to optimize for modern hardware. A key kernel is Sparse General Matrix-Matrix Multiplication (SpGEMM), which underpins simulations, graph analytics, and machine learning applications. SpGEMM exhibits irregular memory access patterns and workload imbalance, making it challenging to achieve high performance on GPUs. Current GPU SpGEMM solutions typically rely on a two-pass workflow to address load imbalance and reduce memory access. The symbolic pass, which determines the number of output elements per row, accounts for roughly 28% of the total runtime on average. In this work, we question the necessity of exact symbolic computation and introduce an estimation-based SpGEMM workflow. Our approach replaces the costly symbolic step with lightweight HyperLogLog estimators, combined with an analysis strategy that dynamically selects the optimal workflow and guides accumulator configuration. In addition, we introduce a hybrid accumulator design, including a novel hash-based accumulator that leverages both shared and global memory. Our approach consistently outperforms leading GPU SpGEMM implementations across a wide range of both square and rectangular matrices, achieving speedups of 1.4x-2.8x on NVIDIA A100 and H100 architectures.
0
0
quant-ph 2026-04-21

Interface abstracts block-encodings for quantum algorithms

by Matic Petrič, René Zander

Block-encodings as programming abstractions: The Eclipse Qrisp BlockEncoding Interface

Constructors and composition tools let users apply these embeddings to matrix inversion and simulation without low-level circuit work.

Figure from the paper full image
abstract click to expand
Block-encoding is a foundational technique in modern quantum algorithms, enabling the implementation of non-unitary operations by embedding them into larger unitary matrices. While theoretically powerful and essential for advanced protocols like Quantum Singular Value Transformation (QSVT) and Quantum Signal Processing (QSP), the generation of compilable implementations of block-encodings poses a formidable challenge. This work presents the BlockEncoding interface within the Eclipse Qrisp framework, establishing block-encodings as a high-level programming abstraction accessible to a broad scientific audience. Serving as both a technical framework introduction and a hands-on tutorial, this paper explicitly details key underlying concepts abstracted away by the interface, such as block-encoding construction and qubitization, and their practical integration into methods like the Childs-Kothari-Somma (CKS) algorithm. We outline the interface's software architecture, encompassing constructors, core utilities, arithmetic composition, and algorithmic applications such as matrix inversion, polynomial filtering, and Hamiltonian simulation. Through code examples, we demonstrate how this interface simplifies both the practical realization of advanced quantum algorithms and their associated resource estimation.
1 0
0
cs.MS 2026-04-09

FAMLIES vertically integrates BLIS and libflame for unified HPC linear algebra

by Devin A. Matthews, Tze Meng Low +3 more

A Proposed Framework for Advanced (Multi)Linear Infrastructure in Engineering and Science (FAMLIES)

The framework targets high-performance linear and tensor operations from single nodes to massive parallelism on CPU and GPU hardware.

Figure from the paper full image
abstract click to expand
We leverage highly successful prior projects sponsored by multiple NSF grants and gifts from industry: the BLAS-like Library Instantiation Software (BLIS) and the libflame efforts to lay the foundation for a new flexible framework by vertically integrating the dense linear and multi-linear (tensor) software stacks that are important to modern computing. This vertical integration will enable high-performance computations from node-level to massively-parallel, and across both CPU and GPU architectures. The effort builds on decades of experience by the research team turning fundamental research on the systematic derivation of algorithms (the NSF-sponsored FLAME project) into practical software for this domain, targeting single and multi-core (BLIS, TBLIS, and libflame), GPU-accelerated (SuperMatrix), and massively parallel (PLAPACK, Elemental, and ROTE) compute environments. This project will implement key linear algebra and tensor operations which highlight the flexibility and effectiveness of the new framework, and set the stage for further work in broadening functionality and integration into diverse scientific and machine learning software.
0
0
cs.MS 2026-04-09 Recognition

Benchmark turns k-server conjecture into inequality search for agents

by Kirill Brilliantov, Etienne Bamas +1 more

k-server-bench: Automating Potential Discovery for the k-Server Conjecture

Agents fully solve k=3 and cut violations on open k=4, giving researchers a new testbed and discovery tool

Figure from the paper full image
abstract click to expand
We introduce a code-based challenge for automated, open-ended mathematical discovery based on the $k$-server conjecture, a central open problem in competitive analysis. The task is to discover a potential function satisfying a large graph-structured system of simple linear inequalities. The resulting evaluation procedure is sound but incomplete: any violated inequality definitively refutes a candidate, whereas satisfying all inequalities does not by itself constitute a proof of the corresponding conjecture's special case. Nevertheless, a candidate that passes all constraints would be strong evidence toward a valid proof and, to the best of our knowledge, no currently known potential achieves this under our formulation in the open $k=4$ circle case. As such, a successful candidate would already be an interesting contribution to the $k$-server conjecture, and could become a substantial theoretical result when paired with a full proof. Experiments on the resolved $k=3$ regime show that current agentic methods can solve nontrivial instances, and in the open $k=4$ regime they reduce the number of violations relative to existing potentials without fully resolving the task. Taken together, these results suggest that the task is challenging but plausibly within reach of current methods. Beyond its relevance to the $k$-server community, where the developed tooling enables researchers to test new hypotheses and potentially improve on the current record, the task also serves as a useful \emph{benchmark} for developing code-based discovery agents. In particular, our $k=3$ results show that it mitigates important limitations of existing open-ended code-based benchmarks, including early saturation and the weak separation between naive random baselines and more sophisticated methods.
0
0
cs.MS 2026-04-09 Recognition

Benchmarks guide MATLAB backend choice for multivariate polynomials

by Yi-Shuai Niu, Shing-Tung Yau

Polylab: A MATLAB Toolbox for Multivariate Polynomial Modeling

MPOLY for light interactive work; MPOLY-HP for heavy simplification and medium-to-large affine-normal tasks; stochastic variant for large-sp

Figure from the paper full image
abstract click to expand
Polylab is a MATLAB toolbox for multivariate polynomial scalars and polynomial matrices with a unified symbolic-numeric interface across CPU and GPU-oriented backends. The software exposes three aligned classes: MPOLY for CPU execution, MPOLY_GPU as a legacy GPU baseline, and MPOLY_HP as an improved GPU-oriented implementation. Across these backends, Polylab supports polynomial construction, algebraic manipulation, simplification, matrix operations, differentiation, Jacobian and Hessian construction, LaTeX export, CPU-side LaTeX reconstruction, backend conversion, and interoperability with YALMIP and SOSTOOLS. Versions 3.0 and 3.1 add two practically important extensions: explicit variable identity and naming for safe mixed-variable expression handling, and affine-normal direction computation via automatic differentiation, MF-logDet-Exact, and MF-logDet-Stochastic. The toolbox has already been used successfully in prior research applications, and Polylab Version 3.1 adds a new geometry-oriented computational layer on top of a mature polynomial modeling core. This article documents the architecture and user-facing interface of the software, organizes its functionality by workflow, presents representative MATLAB sessions with actual outputs, and reports reproducible benchmarks. The results show that MPOLY is the right default for lightweight interactive workloads, whereas MPOLY-HP becomes advantageous for reduction-heavy simplification and medium-to-large affine-normal computation; the stochastic log-determinant variant becomes attractive in larger sparse regimes under approximation-oriented parameter choices.
0
0
cs.MS 2026-04-07 Recognition

Residue overrides eliminate most false reports in floating-point debuggers

by Yumeng He, Pavel Panchekha

Accurate Residues for Floating-Point Debugging

Splitting error computation from residue evaluation and reassembling results across a few runs cuts spurious warnings on scientific code

Figure from the paper full image
abstract click to expand
Floating-point arithmetic is error-prone and unintuitive. Floating-point debuggers instrument programs to monitor floating-point arithmetic at run time and flag numerical issues. They estimate residues, i.e., the difference between actual floating-point and ideal real values, for every floating-point value in the program. Prior work explores various approaches for computing these residues accurately and efficiently. Unfortunately, the most efficient methods, based on "error-free transformations", have a high rate of false reports, while the most accurate methods, based on high-precision arithmetic, are very slow. This paper builds on error-free-transformations-based approaches and aims to improve their accuracy while preserving efficiency. To more accurately compute residues, this paper divides residue computation into two steps (rounding error computation and residue function evaluation) and shows how to perform each step accurately via careful improvements to the current state of the art. We evaluate on 44 large scientific computing workloads, focusing on the 14 benchmarks where prior tools produce false reports: our approach eliminates false reports on 10 benchmarks and substantially reduces them on the remaining 3 benchmarks. Moreover, complex numerical issues require additional care due to absorption, where two machine-precision residues cannot both be computed accurately in a single execution. This paper introduces residue override, which re-executes the program multiple times, computing different residues in different executions and assembling a final "patchwork" execution. We evaluate on 169 standard benchmarks drawn from numerical analysis papers and textbooks, requiring only 3.6 re-executions on average. Among 34 benchmarks with false reports in the initial run, residue override is triggered on 29 of them and reduces false reports on 25 of them, averaging 7.1 re-executions.
0
0
eess.SY 2026-04-03 2 theorems

Hybrid solver reaches exact AVI solutions in finite steps

by Daniel Arnström, Emilio Benenati +1 more

texttt{DR-DAQP}: An Hybrid Operator Splitting and Active-Set Solver for Affine Variational Inequalities

Douglas-Rachford splitting with active-set Newton corrections beats PATH by up to two orders of magnitude on random problems

Figure from the paper full image
abstract click to expand
We present \texttt{DR-DAQP}, an open-source solver for strongly monotone affine variational inequaliries that combines Douglas-Rachford operator splitting with an active-set acceleration strategy. The key idea is to estimate the active set along the iterations to attempt a Newton-type correction. This step yields the exact AVI solution when the active set is correctly estimated, thus overcoming the asymptotic convergence limitation inherent in first-order methods. Moreover, we exploit warm-starting and pre-factorization of relevant matrices to further accelerate evaluation of the algorithm iterations. We prove convergence and establish conditions under which the algorithm terminates in finite time with the exact solution. Numerical experiments on randomly generated AVIs show that \texttt{DR-DAQP} is up to two orders of magnitude faster than the state-of-the-art solver \texttt{PATH}. On a game-theoretic MPC benchmark, \texttt{DR-DAQP} achieves solve times several orders of magnitude below those of the mixed-integer solver \texttt{NashOpt}. A high-performing C implementation is available at \textt{https://github.com/darnstrom/daqp}, with easily-accessible interfaces to Julia, MATLAB, and Python.
1 0
0
cs.MS 2026-03-17 2 theorems

Branch-free algorithms speed up high-precision arithmetic

by Tomonori Kouya

Acceleration of multi-component multiple-precision arithmetic with branch-free algorithms and SIMD vectorization

Multi-component computations using hardware binary64 and binary32 gain acceleration for triple and quadruple precision on x86 and ARM

Figure from the paper full image
abstract click to expand
Multiple-precision floating-point branch-free algorithms can significantly accelerate multi-component arithmetic implemented by combining hardware-based binary64 and binary32, particularly for triple- and quadruple-precision computations. In this study, we achieved benchmark results on x86 and ARM CPU platforms to quantify the accelerations achieved in linear computations and polynomial evaluation by integrating these algorithms.
0
0
cs.DC 2026-03-11 Recognition

FP64 tensor cores speed finite-element kernels 2x

by Jiqun Tu, Ian Karlin +5 more

Accelerating High-Order Finite Element Simulations at Extreme Scale with FP64 Tensor Cores

Kernel fusion on Grace Hopper and Blackwell GPUs yields the gains at 10,000-GPU scale for production codes.

Figure from the paper full image
abstract click to expand
Finite element simulations play a critical role in a wide range of applications, from automotive design to tsunami modeling and computational electromagnetics. Performing these simulations efficiently at the high resolutions needed for practical applications and scientific insights necessitates the use of high-order methods and large-scale supercomputing. While much progress has been made in porting finite element codes to GPU systems in recent years, additional improvements in the efficiency and computational speed of GPU-accelerated high-order finite element simulations are in constant demand. In this paper, we demonstrate that the FP64 tensor cores on NVIDIA GPUs can be used to further accelerate such simulations, achieving significant speedups in key kernels of MFEM, a scalable open-source finite element library widely used in HPC applications. By integrating FP64 tensor cores with kernel fusion optimizations, we were able to achieve up to 2$\times$ performance gains and up to 83% energy efficiency gains on NVIDIA's Grace Hopper GH200 and Grace Blackwell GB200 architectures. To the best of our knowledge, this is the first time that FP64 tensor cores have been directly programmed to accelerate large-scale finite element scientific computing applications. We demonstrate the performance of the optimized kernels at exascale by showing near-perfect weak scaling efficiency and 90% strong scaling efficiency across nearly 10,000 GPUs on the Alps system. The new algorithms and MFEM enhancements directly benefit complex production codes, including the 2025 Gordon Bell Prize-winning application for real-time tsunami forecasting.
1 0

browse all of cs.MS → full archive · search · sub-categories